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Kurzfassung 


Diese Arbeit leistet einen Beitrag zur Optimierung von lokalen endlosfaser- 
verstärkten Patches, unter Berücksichtigung von Fertigungsbedingungen. 


Die Kombination von diskontinuierlich und kontinuierlich faserverstarkten 
Kunststoffen bietet ein breites Anwendungsspektrum. Dabei kann die hohe 
Gestaltungsfreiheit der diskontinuierlich faserverstärkten Kunststoffe 
(DiCoFRP) mit den hohen spezifischen Materialeigenschaften der kontinuier- 
lich faserverstärkten Kunststoffe (CoFRP) kombiniert werden. Dieser Ansatz 
erfordert allerdings spezifische Optimierungsstrategien. Daher wurde eine 
Mehrziel-Optimierungsstrategie, in welcher Fertigungsrandbedingungen Be- 
rücksichtigung finden, für die Positionierung und Dimensionierung von loka- 
len Endlosfaserverstärkungen entwickelt. 


Der in dieser Arbeit entwickelte Mehrziel-Optimierungsansatz verwendet da- 
bei einen evolutionären Algorithmus als Optimierungsalgorithmus. Da diese 
Klasse von Algorithmen eine Vielzahl von Funktionsauswertungen erfordert, 
sind effiziente Methoden zur Bewertung der Fitnesswerte notwendig. Aus 
diesem Grund wird eine kinematische Drapiersimulation zur Vorhersage der 
Patch-Geometrie verwendet. Um die Fähigkeit der kinematischen Drapierung 
zu demonstrieren, wird ein Vergleich mit einem mechanischen Ansatz durch- 
geführt. Dieser Vergleich zeigt den Vorteil der kinematischen Drapiersimula- 
tion, die geringe Rechenzeit, sowie deren Nachteil, eine weniger genaue Vor- 
hersage des Umformverhaltens. 


Als Zielfunktionen für die Mehrzieloptimierung werden sowohl strukturelle 
als auch prozessbezogene Ziele verwendet. Als strukturelles Optimierungsziel 
wird die Bauteilsteifigkeit verwendet, welche mittels einer linear-elastischen 
Struktursimulation ermittelt wird, während für das prozessbezogene Ziel eine 
Verzugssimulation durchgeführt wird. Um die prozessbezogenen Zielfunktio- 
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nen, Bauteilverzug und Eigenspannungen, zu ermitteln, muss das Aushärte- 
verhalten modelliert werden. Hierfür wird eine Abaqus Subroutine vom Typ 
UEXPAN verwendet, mit welcher sich das Aushärteverhalten und die daraus 
resultierenden Dehnungen effizient bestimmen lassen. Neben dem Verzugs- 
ziel werden weitere relevante Fertigungseinflüsse und -randbedingungen aus 
dem Halbzeug-, Drapier- und Co-Molding-Prozess diskutiert und finden wäh- 
rend der Optimierung Berücksichtigung. Die Visualisierung der Optimierungs- 
ergebnisse erfolgt mit Hilfe von Heat-Maps, in denen Bereiche hervorgeho- 
ben werden, welche den größten Einfluss auf die Optimierungsziele haben. 
Darüber hinaus werden Ungenauigkeiten aus dem Herstellungsprozess mit- 
tels einer Robustheitsbewertung berücksichtigt. Hierfür wurde ein Workflow 
entwickelt, um die beiden Robustheitsbewertungsgrößen, Degree of Robust- 
ness und Robustness Index, im Nachgang zur Optimierung zu berechnen. Die 
Berechnung erfolgt mit einem rechnerisch effizienten Ansatz, unter Verwen- 
dung eines Kriging-Metamodells welches auf den bei der Optimierung gewon- 
nenen Daten aufbaut. 


Die entwickelte Optimierungsstrategie wird anhand einer Reihe von Anwen- 
dungsbeispielen, angefangen bei einfachen 2D-Geometrien bis hin zu einem 
komplexen 3D-Beispiel, demonstriert. Dabei wird der Einfluss verschiedener 
Einstellungen des Optimierungsalgorithmus diskutiert. Dabei wird auch der 
Einfluss der Anzahl der Zielfunktionen auf die Performance des Optimierungs- 
ansatzes anhand der 3D-Beispielstruktur demonstriert. 


Abstract 


In this work, contributes to the optimization of local continuous fiber rein- 
forcement patches, under consideration of manufacturing constraints. 


The combination of discontinuous and continuous fiber reinforced plastics of- 
fers a wide range of applications. Thereby the high degree of design freedom 
of the discontinuous fiber reinforced plastics (DiCoFRP) can be combined with 
the high specific material properties of continuous fiber reinforced plastics 
(CoFRP). This approach requires specific optimization strategies. Therefore, 
an multi-objective optimization strategy for the placement of local reinforce- 
ment patches, under consideration of manufacturing constraints, has been 
developed. 


The multi-objective optimization approach developed in this thesis uses an 
evolutionary algorithm as optimization algorithm. Since this class of algo- 
rithms requires a large number of function evaluations, efficient methods for 
the evaluation of the fitness values are necessary. Therefore, a kinematic 
draping simulation is used for the prediction of the final patch geometry. To 
demonstrate the ability of the kinematic draping approach, a comparison 
with a mechanical approach is performed. This comparison shows the ad- 
vantage of the kinematic draping simulation, low computational time, as well 
as its disadvantage, a less precise prediction of the forming behavior. 


During the multi objective optimization, structural and process related objec- 
tives are considered. As structural objective the components stiffness is taken 
from a linear-elastic simulation, while for the process related objective a 
warpage simulation has to be performed. For the prediction of the warpage 
related objectives, distortion and residual stress, the curing behavior has to 
be modeled. An Abaqus subroutine UEXPAN is used for the efficient modeling 
of the curing behavior and the resulting strains. Besides the warpage objec- 
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tive, further relevant manufacturing constraints, coming from the semi-fin- 
ished product, draping process, and co-molding process, are discussed and 
considered as constraints during the optimization. The visualization of the op- 
timization results is done with heat-maps, where areas are highlighted in 
which a reinforcement has the largest impact on the objectives. Furthermore, 
inaccuracies from the manufacturing process are taken into account during a 
robustness evaluation. Therefore, a workflow has been developed to calcu- 
late the two robustness evaluation measures degree of robustness and ro- 
bustness index. The calculation of these measures is done with a computa- 
tional efficient approach, by using a Kriging meta model, which is built on the 
data gained during the optimization. 


The developed optimization strategy is demonstrated with a number of ap- 
plication examples, starting with simple 2D geometries up to a complex 3D 
example. Thereby, the influence of different settings of the optimization al- 
gorithm are discussed. The effect of increasing the number of objective func- 
tions on the optimization performance is demonstrated with the 3D example 
structure. 
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1 Introduction 


1.1. Motivation 


Optimization tools are often part of the product development process, where 
problem specific strategies are necessary to find the best solution. Thereby, 
the range of optimization strategies varies depending on the tools used for 
the initial geometry concept creation, such as topology optimization, to meth- 
ods used in a later steps of the design process, such as shape optimization. 
When using composite materials, the number of design variables is signifi- 
cantly increased, since the material composition itself can serve as an optimi- 
zation objective. In addition, the high anisotropy of this material class has to 
be considered during the optimization process, to fully exploit the excellent 
weight-specific stiffness of composite materials. For this purpose, optimiza- 
tion strategies specifically developed for composite materials are necessary. 


Furthermore, the part’s material cost can be reduced when combining dis- 
continuous fiber reinforced polymers (DiCoFRP) with local reinforcements, 
such as continuous fiber reinforced patches (CoFRP), to create load carrying 
paths. Therefore, acommonly used optimization objective is to achieve high- 
est performance in terms of maximum stiffness, while using a minimum 
amount of reinforcement material. In addition, the concurrent use of contin- 
uous and discontinuous fiber reinforced composites, allows a combination of 
the advantages of both material classes, i.e. DiCoFRP and CoFRP. While the 
excellent weight specific properties of the continuous reinforced material are 
utilized in highly loaded areas, the discontinuous reinforced composite mate- 
rial offers a wide range of design freedom. To achieve an advantageous com- 
bination of both material classes, an optimization strategy dealing with the 
positioning of the continuous fiber reinforcements (patches) is necessary. 
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In general, the use of virtual optimization tools reduces the number of expen- 
sive prototypes, since the product quality can be increased in an early phase 
of the design process. Therefore, a number of different optimization ap- 
proaches have been developed in the past, mostly focusing on isotropic ma- 
terials. The methods presented for composite optimization are usually focus- 
ing on either continuous or discontinuous reinforcements, while the 
combination of both material classes has not been a research focus. Hence, a 
method should be proposed, to combine both material classes, utilizing the 
advantages of each class. In addition, the quality and feasibility of optimiza- 
tion results should be increased by the incorporation of manufacturing con- 
straints during the optimization. An additional improvement of the optimiza- 
tions result can be achieved by performing a robustness evaluation. Thereby, 
further influences of the manufacturing process can be incorporated, to en- 
sure that the predicted optimum remains an optimal solution in the real pro- 
cess. 


1.2. Objectives of the thesis 


The global objective of this thesis is the development of a strategy for the 
optimization of position and size of continuous fiber reinforced patches, un- 
der consideration of manufacturing constraints and warpage. 


Therefore, a multi-objective optimization method has to be implemented, 
which allows a parametric optimization of the patches and the determination 
of the final Pareto front, where a Pareto front defines a multi-objective trade- 
off curve. In addition to the Pareto front, a second tool for the designer should 
be proposed to support the patch positioning in the final part. Manufacturing 
constraints, like maximum patch size or feasible patch geometries, shall be 
considered during the optimization. For this purpose, an efficient draping sim- 
ulation has to be implemented in the optimization workflow. Furthermore, a 
warpage objective should be integrated in the optimization workflow to limit 
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arising distortions and inner stresses. Therefore, curing and warpage simula- 
tion models have to be implemented. Draping method and curing simulation 
need to be implemented in a CAE chain to model the entire process from 
starting material to final product, which is finally used for the optimization. In 
addition, a method for the visualization of the optimization results has to be 
developed, which can assist the designer in finding the area, in which the 
patch reinforcement has the largest impact. 


Since the robustness of an optimization result is an important objective, a 
method should be proposed, which allows for evaluation of the final front’s 
robustness. Therefore, a meta-model based robustness evaluation method 
has to be implemented, which uses the results from the multi-objective opti- 
mization as input information. 


1.3. Structure of the thesis 


The thesis is organized as follows: An overview of the state of the art, regard- 
ing the methods used and developed in this thesis is given in Chapter 1. In 
Chapter 2, the patch optimization method is introduced, and an introduction 
in the field of evolutionary algorithm is given, including basic definitions re- 
garding optimization. 


The used kinematic draping simulation method is presented in Chapter 3. In 
addition, a comparison of kinematic and mechanical draping methods is given 
here. The approach of the proposed curing simulation is introduced in Chap- 
ter 4. Besides the curing method, the approach used for the calculation of the 
residual strains, occurring during cure, is given here. 


The optimization along the CAE chain, developed within this work, is pre- 
sented in Chapter 5. In addition, the constraints, resulting from the manufac- 
turing process, to be considered within the optimization are discussed in this 
chapter. The developed optimization approach is applied to a number of ap- 
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plication examples, numerically verified and discussed in Chapter 6. To fur- 
ther include manufacturing inaccuracies in the optimization process, a robust- 
ness evaluation method is developed. The proposed workflow of the robust- 
ness evaluation method is presented in Chapter 7. First, an overview of 
robustness evaluation, in the context of multi-objective optimization is given. 
Furthermore, a surrogate model approach, used for the robustness evalua- 
tion, is presented. In addition, the developed method is applied to an appli- 
cation example, introduced in the previous chapter. 


This thesis is concluded with Chapter 8, summarizing the findings and giving 
an outlook for further research in the field of composite optimization. 


2 State of the Art 


2.1. Overview of the SMC process with continuous 
fiber reinforcements 


The optimization strategy of this thesis focuses on the compression molding 
process of sheet molding compounds (SMC) with over-molding of continuous 
fiber reinforced patches. Therefore, a brief introduction in the processing of 
patch-reinforced SMC parts is given here. 


The first step is the preparation of pre-impregnated semi-finished products 
(prepreg) with discontinuous and continuous fiber reinforcement (DiCoFRP 
and CoFRP), cf. Figure 2-1 (a) and (b). After the resin is applied to the bottom 
foil, either continuous or discontinuous fibers are added, and finally a resin- 
coated top foil is applied. Within the GRK2078 project frame of this thesis, an 
unsaturated polyester polyurethane hybrid (UPPH) resin system is used for 
both prepregs [1]. Once the fabrication of the semi-finished products is done, 
the materials are cut and stacked. After a near net-shaped draping process of 
the CoFRP material, cf. Figure 2-1 (c), the material is positioned in the press 
and co-molded, cf. Figure 2-1 (e). The last steps of the process are the demold- 
ing and deburring of the final parts, cf. Figure 2-1 (f). The reader is referred to 
Sharma et al. [2] for further details on the SMC process, and to Bücheler et al. 
[1] for the specific requirements of the combination of SMC with local patch 
reinforcements. The manufacturing of components with multiple fiber archi- 
tectures, like SMC with unidirectional reinforcements, and their behavior dur- 
ing the manufacturing process, was also discussed by Corbridge et al. [3]. 


The lightweight ability of the combination of continuous and discontinuous 
fiber reinforced materials has been demonstrated within several industrial 
applications. Akiyama [4] presented a trunk lid and an engine cover, while 
Bruderick et al. [5] showed a more holistic approach with the design of the 
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Dodge Viper, where several SMC components are used in combination with 
CoFRP. 


Figure 2-1: Process for continuous-discontinuous-fiber reinforced thermoset components, 
used within GRK2078; discontinuous-fiber prepreg preparation (a), continuous fiber 
prepreg preparation (b), handling and near net-shape draping (c), position and ori- 
entation control of the continuous reinforcement (d), co-molding (e), and final de- 
burring (f) [1] 


2.2. Structural optimization 


In product development, optimization methods have gained an important 
role, since they allow a designer to improve a component’s performance in 
early design phases. Therefore, several different approaches have been de- 
veloped in the past. Structural optimization methods can be classified in the 
three main groups sizing, shape, and topology optimization. For these ap- 
proaches the degree of design freedom increases from sizing to topology op- 
timization. While using a size optimization method, only a given set of dimen- 
sions are subject to the optimization, such as the cross section thickness of 
beam elements, cf. Figure 2-2. Whereas, shape optimization is a more flexible 
approach, allowing also to change the shape of the geometry. Here usually a 
set of control points for morphing is defined to describe a flexible geometry. 
These points work then as degrees of freedom for the optimization. In the 
example given in Figure 2-2, the shape of the holes would be subject to the 
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optimization process. The approach with the highest degree of freedom, re- 
garding geometry changes, is the topology optimization. In a product devel- 
opment process, topology optimization is usually used in an early state to cre- 
ate a conceptual design [6]. While size and shape optimization start from an 
almost final part design, topology optimization starts from a rough design 
space with a uniform material distribution. Therefore, usually a volume con- 
straint is used, to describe the maximum allowed amount of material [7]. For 
a detailed overview on topology optimization methods, the reader is referred 
to Eschenauer and Olhoff [8]. A comprehensive overview and evaluation of 
structural optimization methods in general was published by Saitou et al. [6]. 
The approach presented here is a combination of shape and topology 
optimization, where the change of patch size is a kind of shape optimization, 
while the positioning of the patch can be assigned to the field of topology 
optimization. 


je 
Sizing optimization 


CEE KK) - Cx) 


Shape optimization 


a a LLNS 
Topology optimization 


Figure 2-2: Comparison of structural optimization approaches [9] 


Since for most applications the quality of optimization results can be in- 
creased by the incorporation of manufacturing constraints, this topic has 
been a field of research for some time. One of the first constraints to be in- 
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cluded in topology optimization was the consideration of an extrusion direc- 
tion, which ensures that the resulting geometry can be manufactured by an 
extrusion process [10]. In the same work, the authors developed an approach 
for the consideration of a mold opening direction in the topology optimiza- 
tion. While using a cast molding part as an application example, Zhou et al. 
demonstrated the influence of manufacturing constraints on the quality of 
the optimization result. 


Zuo et al. [11] introduced further machining and manufacturing constraints 
for the topology optimization. To ensure a post-processing, such as milling, 
they proposed a minimum hole size control, which allows them to remove 
small holes form the resulting structure. Furthermore, they integrated a sym- 
metry constraint in the optimization, which is able to ensure demoldability 
for injection molded parts. 


The approaches from Zhou et al. [10] and Zuo et al. [11] were integrated in 
the so-called Solid Material with Penalization Model (SIMP), where a virtual 
density is used as optimization parameter. Therefore, in both approaches the 
density distribution over the thickness of the model is used as constraint, in a 
way that no change of the gradient in thickness direction is allowed. For a 
detailed explanation of the SIMP approach the reader is referred to Harz- 
heim [7]. 


The importance of manufacturing constraints for industrial applications is 
demonstrated by Gruber et al. [12]. They used topology optimization with in- 
tegrated minimum member size control, to generate design proposals for an 
aircraft structure. Liu et al. [13] discussed the influence of different manufac- 
turing constraints on the topology optimization results. For this purpose, they 
compared different constraints like minimum hole size (to prevent small 
holes) or length control methods (to prevent checker boarding). Focusing on 
injection molded parts, they demonstrated the capability of a rib thickness 
gradient method to ensure the same thickness along a rib, which is useful to 
prevent problems with warpage. They also demonstrate the usefulness of the 
extrusion and tool opening direction constraint. 


2.2. Structural optimization 


Usually gradient-based optimization approaches are used for topology opti- 
mization. Chapman et al. [14] showed in their work the integration of an evo- 
lutionary algorithm in the topology optimization. In their approach, they use 
the elements as design variables, which results in a large number of design 
variables, even for small models. Thereby, every element corresponds to an 
entry in the vector of design variables, which is in the context of evolutionary 
algorithm also called gen string. For larger structures, they proposed an ap- 
proach for the building of sub-groups as well as a connectivity check. Within 
their work they focused only on 2D structures. 


Even though size and shape optimization are usually used in later steps of the 
design process, those methods can also be used to provide input for initial 
design concepts. For this purpose a truss structure, consisting of beam ele- 
ments, is used as a representation for the geometrical optimization problem. 
Those representations are typically used in civil engineering optimization 
problems. [15] 


Like topology optimization, shape optimization has been a topic in research 
for a while. Haftka and Grandhi [16] give a good introduction of the applied 
methods. Furthermore, they describe some reccurring problems, like stress 
localization due to sharp edges generated by the change of shape or the need 
of mesh refinement during the optimization. An overview of the latest devel- 
opments and trends, with focus on aerodynamic applications, for shape opti- 
mization methods was given by Skinner and Zare-Behtash [17]. It has been 
proven by several authors, that including manufacturing constraints in the 
optimization workflow significantly increases the quality of the results [11- 
13]. In addition, it could be seen that evolutionary algorithms are able to deal 
with complex optimization problems, like topology optimization [14]. 


Beyond the described methods, approaches with the integration of surro- 
gate-models or reduced-order models have been used for structural optimi- 
zation [6]. Additionally an approach to use particle swarm optimization (PSO) 
methods as algorithm for a structural optimization problem have been pre- 
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sented by Perez and Behdinan [18]. A PSO approach for 2D topology optimi- 
zation has been presented by Luh et al. [19], thereby they used the elements 
as discreet design parameters. Another field gaining attention for the devel- 
opment of optimization methods is the additive manufacturing, since this 
technology is able to realize relatively complex structures, following mostly 
the result from the optimization. Therefore, Liu et al. [20] addressed this topic 
in their work, by the consideration of deposition paths during the optimiza- 
tion run. Contrary to evolutionary algorithm, PSO is usually used for single 
objective optimization only, furthermore, a rather complex modeling of the 
optimization problem is necessary, as input for the optimization process. 


2.3. Composite optimization 


Since composite materials offer several additional degrees of design freedom, 
such as fiber or matrix type, fiber volume content, or stacking sequence of 
laminates, special optimization strategies are necessary to fully exploit the 
abilities of this material class. Common topology optimization approaches are 
not able to cover all possible degrees of design freedom. Therefore, several 
approaches have been made to expand the capability of topology optimiza- 
tion for composite applications. 


Topology and structural optimization for composite applications 


In their work, Blasques and Stolpe [21] utilise the classical laminate theory 
(CLT) within a topology optimization. Based on a predefined set of layer types, 
a list of possible layups, each consisting of several layers, is created a-priori. 
Those materials where then assigned during the topology optimization, 
thereby the material assignment represents an additional parameter, com- 
pared with the common topology optimization approaches. Dai et al. [22] ex- 
tended this approach by integrating the CLT directly in the optimization loop, 
considering 2D geometries only, but offering more degrees of freedom for 
the optimization. Therefore, not only the predefined layups can be used, but 
all possible layer combinations. 
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A more holistic approach is presented by Liu et al. [23], where they use three 
different objectives, resulting from a static load case, a frequency analysis, 
and a crash simulation. Therefore, the model is divided by a user-predefined 
segmentation, each of this segments has a specific layup and thickness, where 
the thickness serves as design parameter. For these parameters, a specified 
number of simulations is done to create a Kriging surrogate model, which is 
then used as input for a particle swarm optimization algorithm. 


As the fiber orientation has a major influence on the mechanical properties, 
orientation optimization strategies have been addressed by several authors 
[24-28]. Savic et al. [27] used an improved hit and run algorithm for the opti- 
mization of fiber orientations in I-profiles. Therefore, they subdivided the pro- 
file in several regions and used the local fiber orientation of each region as a 
design parameter. Focusing on 2D structures, Legrand et al. [28] proposed a 
method for the optimization of the element orientation. Here, they inte- 
grated a genetic algorithm in a finite element simulation, with the global 
strain energy as objective function. Extending this approach, Vosoughi et al. 
[25] used a genetic algorithm combined with a particle swarm method to 
maximize the buckling load of a given structure. 


Process parameter optimization 


Since the mechanical properties of composite materials are strongly depend- 
ent on the manufacturing process, it is self-evident to study the abilities of 
virtual process optimization. An approach for the minimization of the shear 
angle is presented by Skordos et al. [29], thereby they used a combination of 
an evolutionary algorithm and a kinematic draping simulation, with the ori- 
entation of the layers as an optimization parameter. In addition, an optimiza- 
tion of the draping process with a mechanical-based draping simulation was 
proposed by Chen et al. [30], where they used a series of springs to model 
and optimize the gripping system used during the draping process. A two-loop 
CoFRP optimization workflow, considering a process and structural simulation 
has been developed by Karger et al. [31]. In their work, they used the process 
simulation to optimize the draping process, while the structural performance 
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is evaluated subsequently by a structural simulation, considering the results 
from the process optimization. A further method to consider draping effects 
in the optimization is introduced by Kaufmann et al. [32]. Since they focus on 
the weight optimization, they create an a-priori draping data base, which then 
can be used during the optimization loop. This means the draping simulation 
is not part of the optimization loop. 


With the minimization of residual stresses, Sonmez and Akbulut [33] focus on 
the optimization of the tape laying process itself. Therefore, they included a 
thermoplastic tape laying process simulation in the optimization loop, while 
searching the optimum process parameters. In the optimization workflow, 
the tape laying process simulation is combined with a residual stress analysis, 
including a degradation model, and a simulation of the bonding process. 
Therefore, they used the process parameters such as temperature, roller ve- 
locity, and heated length as design variables. With this, they performed a sin- 
gle objective optimization to identify the optimal process window, in terms of 
residual stress minimization, with the degree of bonding and the thermal deg- 
radation as constraints. 


Optimization of layup and stacking sequence 


A topic covered by several authors is the search for the optimum layup con- 
figuration. This has been first addressed by Le Riche and Haftka [34], where 
they used an evolutionary algorithm to optimize the buckling load of a struc- 
ture. An extended approach was introduced by Zehnder and Ermanni [35], 
with their patch based optimization of a boat structure. In their context, a 
patch is defined as an area with a constant layup and thickness, within a struc- 
ture. These patches have to be defined a priori, and their shape will remain 
the same during the optimization, while the stacking sequence in each patch 
is subject to the optimization. 


The latest development has been presented by Albanesi et al. [36], where 
they extended the optimization quality by including manufacturing con- 
straints, besides mechanical constraints, to minimize a component’s weight. 
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Therefore, they combined a genetic algorithm with a Laplace crossover and a 
power mutation approach, using ply order, ply number, and ply drop as de- 
sign parameter. They used a wind turbine blade as reference model, to 
demonstrate the proposed approach, cf. Figure 2-3. 


EB] gelcoat 
12 E] uniaxial fiberglass 
a - biaxial fiberglass 
9 n T 7 EEE double-bias fiberglass 


Plies 


Patches 


Figure 2-3: Example of a stacking sequence optimization for a wind turbine blade [36] 


Since composite materials consist of several components, there has been 
some work done, dealing with the optimization of the material composition. 
Sadagopan and Pitchumani [37] suggested a method to minimize manufac- 
turing cost, while maintaining mechanical and thermal properties. In their 
work, they compared a genetic algorithm and a simulated annealing ap- 
proach, while searching for the optimum material composition, in terms of 
fiber architecture, fiber and matrix type, and fiber volume fraction. 


A special case of reinforcing structures is covered by Jansson et al. [38] with 
the optimization of local unidirectional reinforcement tows, cf. Figure 2-4, 
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where the authors proposed a method of the combined optimization of the 
layup and tow thickness, while the geometry remains unchanged. 


Figure 2-4: Example of a local tow reinforcement [38] 


Rettenwander et al. [26] introduced an iterative approach for the optimiza- 
tion of local part thickness and fiber orientation. Here, thickness and fiber 
orientation are not modified at the same time, but instead optimized sequen- 
tially after each other. This loop is run through several times. The optimization 
of local thickness and fiber volume content is addressed by Qian et al. [39] 
with an element based approach. To consider manufacturability, they in- 
cluded a combination step, in which they combined small regions with similar 
thickness to larger regions with equal thickness. Thereby, areas which are too 
small to manufacture are removed. Using a spare wheel housing as case 
study, the authors demonstrated the capacity of their approach [40]. 


A method to consider uncertainties during a laminate optimization procedure 
was presented by Kalantari et al. [24], where they used fiber orientation, fiber 
volume fraction, local thickness, and stacking sequence as parameter. In the 
proposed method, they incorporated a lower and an upper bound for the 
consideration of scattering within the three parameter fibers volume fraction, 
fiber orientation, and local thickness. Therefore, in each iteration, an inter- 
mediate step is included, in which a worst-case scenario for the three scatter- 
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ing parameter is assumed, while the total thickness and volume fraction re- 
main constant. The result from this worst-case scenario is then used during a 
stacking sequence optimization. 


While looking at commercially available solutions, the most common method 
is the implementation in Altair’s software OptiStruct, cf. Figure 2-5. Here, a 
procedure consisting of the three main steps free-size optimization, compo- 
site size optimization, and composite shuffling is implemented. The result of 
the free-size step are the required thicknesses of each layer with a defined 
fiber orientation. Based on this result, the Size Optimization is used to calcu- 
late the necessary thickness of each ply, thereby discrete thicknesses are con- 
sidered to ensure manufacturability. In the final step, the stacking sequence 
of the laminate is optimized. 


Figure 2-5: Composite optimization procedure implemented in OptiStruct (Source: 
altairhyperworks.com) 
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Optimization of local patch reinforcements 


In their work, Mathias et al. [41] focused on the optimization of composite 
repair patches for aluminium structures. Therefore, they used a spline-based 
representation of a 2D patch as a reinforcement for damaged structures, cf. 
Figure 2-6. Beside the control points, they used the global patch orientation 
as well as the local fiber orientation as optimization parameter. Zehnder and 
Ermanni [35] proposed a similar approach, by using a parametric CAD model 
to describe the patch geometry for the optimization. 


Figure 2-6: Example of a spline-based definition of a local repair patch [41] 


Rettenwander et al. [42] introduced a method for the optimization of local 
unidirectional reinforcement patches. The positioning of the patches is done 
with respect to a stress field resulting from a finite element simulation, where 
the patches follow the main stress paths. This procedure is based on the CAIO 
method (Computer Aided Internal Optimization), where tailored fiber rovings 
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are oriented along the principle stress trajectories. Therefore, they focused 
on 2D structures and did not include a draping step. Another approach for the 
optimization of the fiber orientation based on the CAIO method was proposed 
by Klein et al. [43]. In their work, they also focused on 2D structures only. 
Focusing on repair strategies, Kashfuddoja and Ramji [44] used a genetic al- 
gorithm for the optimization of the shape of local patches, therefore the di- 
mensions of patches with predefined geometries are used as parameters. The 
optimization is performed with four predefined layups. The design parame- 
ters of the composite patch, like fiber orientation or fiber volume content, 
remain unchanged. 


A different definition for a patch is given by Zhang et al. [45], in their approach 
a patch refers to an area with a uniform stacking sequence. Here, they subdi- 
vided the initial geometry a priori in patches, for which they used the ply ori- 
entation as parameter and a weighted sum, consisting of the two objectives 
compliance and eigenfrequency. 


Optimization of local reinforcements for composite parts has shown its po- 
tential in several different applications. In most cases, manufacturing effects 
are not considered during the optimization process, though these effects can 
largely influence the optimization process, and should therefore be included 
in the workflow. Furthermore, several approaches are only applicable for 2D 
applications, and do not consider draping effects in the structural optimiza- 
tion. In addition, the presented approaches focus mostly on the optimization 
of single objectives, from either a process or structural point of view. Here a 
more holistic approach is necessary, to fully exploit the advantages of com- 
posite materials. 
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2.4. Multi objective optimization 


In general a multi-objective optimization problem can be formulated as 


Minimize [F(x)] = [F; (x), F)(x), Fe)", (1) 
xEX 
subject to 
go j=1,2,.m (2) 
h(x) =0, i=1,2,..n (3) 


where k is the number of objective functions F(x) and x represents the vec- 
tor of design variables in a feasible design space X. Furthermore, the con- 
straints are given by j inequality constraints g(x) and i equality constraints 
h(x). The objective of the optimization is to find a solution x which minimizes 
all objective functions simultaneously. In general, a minimization is sufficient 
for all problems, since a maximization can be formulated as 


Maximize [F(x)] = Minimize [-F(x)]. (4) 


For most engineering applications, the objective functions have opposite cor- 
relations, which means that improving one objective will impair another one. 
To quantify these differing correlations, Villfredo Pareto described the idea to 
find a set of optimal solutions, the so called Pareto set or Pareto front [46- 
48]. The Pareto optimality is defined as a set of solutions x = [x1, X2,..., Xp] 
for which no point x* exists, with F;(x) > F;(x*) for all F;(x) and 
F;(x) > F;(x*) for at least one function F;(x), with i = 2 ...n as the number 
of considered fitness functions. To solve this problem several different ap- 
proaches exist in literature. Thereby, two quality measures are pursued to 
evaluate the optimization algorithms. The first is to ensure convergence in 
the direction of the best solutions, and the second is to maintain a certain 
diversity along the front, to obtain a better representation of the front, cf. 
Figure 2-7. Additionally should be mentioned, that it cannot be ensured that 
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an optimization converges to a global optimum rather than a local. This ap- 
plies in general for multi-objective as for single-objective optimization prob- 
lems. 
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Figure 2-7: Main objectives of a multi-objective optimization: Convergence to an optimal 
front and a high diversity along the final front 


For the most common multi-objective optimization method, the weighted 
sum approach, different strategies are presented in literature [48, 49]. The 
idea of this method is to combine all objective functions on a single objective 
function, in a way that only one function U has to be minimized. For the sim- 
plest case U is defined as 


k 
U= wi Fo), Fix) > Ovi (5) 


i=1 


In which w; represents the individual weighting factors for each objective 
function F;. Several studies have been conducted with these methods, and 
several variations have been published. For the creation of a Pareto set, the 
optimization has to be done with several different sets of w. Furthermore, 
those methods are susceptible to user influence, since the different sets of w 
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have to be selected by the user. Therefore, no non-problem-specific strategy 
exists, which means that the user has to have a good knowledge of the be- 
havior of the objective functions to define a reasonable set of w. [48] 


In addition several other strategies, like the weighted global criterion or the 
lexicographic method, have been developed with the aim to find the Pareto 
optimal set. All of those methods have in common, that multiple optimization 
runs, each with adjustments made by the user, are necessary [48], except for 
genetic algorithm methods. With these methods it is possible to calculate the 
complete Pareto front within one optimization loop. Therefore, these meth- 
ods are insensitive to user influence, since here no weighting vector w has to 
be defined. In the following, an overview of engineering applications of multi- 
objective optimizations is given. Despite this drawback, in most cases a 
weighted sum approach is implemented to solve these problems. 


Walker and Smith [50] used a genetic algorithm for the simultaneous optimi- 
zation of mass and deflection. Therefore, they used a weighted sum method, 
while a Tsai-Wu failure criteria is incorporated as a constraint. Thereby the 
fiber orientation of each ply and the laminate thickness are implemented as 
discreet design parameters. Almeida and Awruch [51] investigated the influ- 
ence of different methods for the representation of the design variables, 
while using a genetic algorithm with a weighted sum. Hemmatian et al. ap- 
plied a gravitational search algorithm [52] and an ant colony optimization 
method [53] to solve a composite optimization problem, with the objectives 
of total weight and cost, while combining different material types. 


A patch optimization procedure is presented by Mejlej et al. [54], where they 
used a combination of two different genetic algorithms. In their work, they 
covered the stiffness and weight optimization, combined in a weighted sum, 
of local reinforced metallic components, but did not include forming or warp- 
age simulation to predict the realized patch position or to consider warpage. 
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Although most references use a weighted sum method, some approaches 
have been done to use Pareto optimization approaches in engineering opti- 
mization. Vo-Duy et al. [55] used the contrary optimization problem, minimi- 
zation of weight and maximization of natural frequency, to demonstrate the 
ability of genetic algorithm in finding the Pareto front. Thereby, they also 
showed the influence of the number of optimization design variables by var- 
ying them from one variable to three variables. Ermakova and Dayyani [56] 
applied a weighted sum method and a genetic algorithm to a shape optimiza- 
tion problem. In their work, they pointed out the problems with the random 
character of genetic algorithm, which in their case results in non-smooth re- 
sults. 


Another field for the application of multi-objective optimization techniques 
are optimization procedures covering several different simulation methods 
for the calculation of the fitness functions. Eck et al. [57] introduced a method 
to consider a process cost objective, resulting from the resin transfer molding 
(RTM) process, within the optimization. Therefore, they used a surrogate 
model tool, which they called process-estimator, to predict the mold filling 
time and total process cost depending on the stacking sequence. A simulta- 
neously structural and aerodynamic optimization approach is introduced by 
Dal Monte et al. [58], where they concurrently optimize the stiffness and the 
generated power of a wind turbine blade. Akmar et al. [59] suggested a multi- 
scale approach for the optimization of hybrid composite layup structures. 
Their approach consists of a genetic algorithm, which is used to optimize the 
weave pattern by using a representative volume element (RVE), and an ant 
colony method for the optimization of the stacking sequence. 


Pohlak et al. [60] developed an approach for the optimization of large com- 
posite parts, where they use cost and time as objectives, while considering 
stress, displacement, and maximum thickness constraints. The optimization 
is performed on a meta-model, which is trained with input data from a free- 
size optimization. 
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It has been proven by several authors, that evolutionary algorithms are able 
to solve complex engineering optimization problems. Nevertheless, most ap- 
proaches utilize a weighted sum approach rather than a multi-objective opti- 
mization, which results in several optimization runs, from which each is sub- 
ject to user influence. Furthermore, it could be seen that the combination of 
different simulation methods is a promising field. Therefore, a structural sim- 
ulation will be coupled in this work with process simulation, warpage predic- 
tion and kinematic draping simulation, to form a multi-objective optimization 
approach. In addition, the use of meta-models has been proved to be useful 
for composite optimization problems. Therefore, a meta-model approach will 
be used for the evaluation of robustness. 


2.5. Draping simulation 


Since the fiber structure determines the mechanical properties of the final 
part, it is important to predict the fiber architecture of continuous fibre rein- 
forcements by draping simulation. In principle, the available draping simula- 
tion approaches can be separated in two fields, kinematic (also known as 
mapping approaches) and mechanical methods [61]. Both approaches have 
been developed with a different focus. Mechanical approaches are based on 
constitutive models and offer higher prediction accuracy of the actual mate- 
rial behavior [62-65]. Therefore, a good knowledge of the material behavior 
is essential. Compared to this, kinematic approaches offer only an approxi- 
mate prediction of the material behavior and, thus, do not require material 
characterization. Moreover, kinematic approaches are usually very efficient. 
Due to the lack of exact modeling of the material behavior, it is not possible 
to predict wrinkling by kinematic draping simulation. Nevertheless, van West 
et al. [66] proposed a knowledge-based prediction of wrinkling by selecting 
the initial paths in a way, that the shear angle can be used as an evidence. 
Therefore, they placed the initial paths close to the point of interest. Never- 
theless, this approach is only feasible for simple geometries, since the initial 
paths have to be placed individually for each possibly problematic geometry 
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feature. Long and Rudd [67] used the kinematic draping simulation to predict 
hotspots of fiber misalignment. Thereby, they used different angles for the 
initial paths as representation of different fiber orientations of the undraped 
fabric. Since the draping simulation used here is a kinematic method, this 
chapter will focus on kinematic draping simulation. 


The basic assumption, made for kinematic draping simulation of unidirec- 
tional fiber reinforcements, is a simple shear deformation of the material, 
where the length of the fibers and the spacing between the fibers is assumed 
to be constant during the forming process, cf. Figure 2-8. The kinematic drap- 
ing approach has been first introduced by Mack and Taylor [68] for the form- 
ing behavior of woven textile fabrics. Woven textiles can be assumed to de- 
form in pure shear, where the length of the fibers in both directions is kept 
constant. The following assumptions are made for kinematic draping simula- 


tions: 
1. fibers are inextensible, 
2. warp and weft crossing points act as pivoting points, 
3. noslipping occurs at the pivoting points, 
4. the distance between two pivoting points is much less than the 


smallest curvature radius of the surface, 
5. the fabric is everywhere in contact with the surface. 


Those constraints serve as basis for every implementation of a kinematic 
draping approach, and will therefore also be applied here. 
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Figure 2-8: Assumptions of the simple shear deformation behavior of unidirectional fiber rein- 
forcements, used for the kinematic draping simulation, where fiber length and dis- 
tance between the fibers remain unchanged during the draping procedure [61] 
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The drawbacks and advantages of kinematic methods have been extensively 
studied by multiple authors. A main disadvantage of kinematic draping simu- 
lation is the dependence of the simulation result on the choice of the starting 
point and of the initial path of the kinematic algorithm. Van West et al. [66] 
demonstrated the influence of the initial path for the draping of a corner 
bending, where the results depend heavily on the selected starting point. 
Vanclooster et al. [69] systematically studied the influence of different orien- 
tations of the initial path and different positions of the starting point. 


Sharma and Sutcliffe [70] used a double curved structure to compare the re- 
sults of a forming experiment with the result from a draping simulation. Fur- 
thermore, they showed that the kinematic draping simulation has problems 
to predict the shear angle in areas with fast changing curvature. Pickett et al. 
[71] used a rear seat bench and wheel arches of an automotive structure to 
prove that a kinematic draping simulation can be used to get a first impression 
of the draping behavior in an early state of product design, since it is reason- 
ably robust and requires only a minimum of data input. Nevertheless, they 
also pointed out that the starting point and initial path has a significant influ- 
ence on the result of the simulation. 


It has been proven by several authors, that a kinematic draping simulation is 
sufficient to give a coarse representation of the patch’s geometry. However, 
a comparison of the results gained with a FE based draping simulation and a 
kinematic draping simulation, is performed within this thesis. Furthermore, 
the kinematic draping simulation has significant benefits, in terms of compu- 
tational time. Since the proposed optimization procedure requires a large 
number of calculations, a kinematic draping approach is used in this work. 


2.6. Curing and warpage simulation 


The causes for warpage during composite manufacturing processes can be 
split up into a thermal and a chemical part. The thermal proportion results 
from the different coefficients of thermal expansion of fiber and matrix, while 
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the source for the chemical proportion is the creation of covalent bonds dur- 
ing the curing of the matrix. The distance between the monomer molecules 
in the liquid state is dominated by Van der Waals interactions. Since the spe- 
cific distance of covalent bonding is about an order of magnitude smaller, the 
curing of the matrix results in a shrinkage of the matrix material. The amount 
of shrinkage is thereby influenced by type and concentration of matrix, hard- 
ener, inhibitor, and fillers [72]. 


The model used in this thesis to describe the curing behavior was first intro- 
duced by Kamal [73], which is a further development of a model proposed 
earlier by Kamal and Sourour [74] including additional fitting parameters. The 
governing equations are presented in Section 5.1. A model for the curing of 
sheet molding compound (SMC) during the compression molding process was 
proposed by Lee [75]. The approach proposed by Lee covers the kinetic mech- 
anism with governing equations for initiation, inhibition and propagation of 
the curing process. Furthermore, Lee presented a method for the simulation 
of the heat transfer within the SMC material, during the compression molding 
process, and compared his simulation results with experimental data. Since 
the chemical composition of the SMC material used in this thesis is not avail- 
able, the method from Lee could not be used here and, thus, the Kamal model 
is chosen. However, the choice of the curing model does not affect the prin- 
cipal goal of the thesis, since the curing simulation is not the main focus, but 
its implementation in the optimization workflow. 


An overview of different kinetic models is given by Halley and Mackay [76], 
where they present the typically used approaches, like first and n* order, pol- 
ynomial, or autocatalytic models. Furthermore, they present the common 
methods for the chemo-rheological characterization, as well as for the cure- 
dependent modeling of the viscosity. 


In their work, Tseng and Osswald [77] presented an FE based method for the 
prediction of shrinkage and warpage of short fiber reinforced thermoset com- 
posites during the compression molding process. Furthermore, they consid- 
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ered the fiber orientation from the process simulation for the shrinkage cal- 
culation. A separation of different effects, leading to warpage, was done by 
Radford and Rennick [78]. According to their work, the most important influ- 
ences are material anisotropy, temperature gradient through parts thickness 
direction, and part-tool interactions. 


Bogetti and Gillespie [79] presented an approach for the prediction of process 
induced stress and deformations within thick-section thermoset composite 
laminates. For this purpose, they proposed a combination of a one-dimen- 
sional cure simulation to calculate the gradients of degree of cure and tem- 
perature over the thickness. A micro-mechanical model is used to calculate 
the temperature and cure dependent material properties, while an incremen- 
tal laminated plate theory is used to calculate the laminate properties, and 
perform the warpage calculation.. In addition, they presented a method to 
calculate the development of the Young’s modulus depending on the degree 
of cure. Based on this work, Svanberg and Holmberg [80-82] developed an 
approach for the calculation of the strains resulting from curing. In their ap- 
proach they split the occurring strains in a chemical and a thermal compo- 
nent. The governing equations are presented in Section 5.2. The presented 
approach is relatively computationally inexpensive, and therefore used within 
this thesis. 


While focusing on unsymmetrical laminates, Wisnom et al. [83] analyzed the 
basic mechanism causing residual stress and distortion, and especially their 
development during cure. Thereby, they showed that after vitrification, the 
main cause for distortion is thermal stress, while before vitrification thermal 
and chemical stresses, have to be considered. In their work, Wisnom et al. 
focused on experiments, and did not perform simulations. Hossain [84] pro- 
posed the inclusion of a modeling of the material’s stiffness increase during 
curing in combination with a phenomenological viscoelastic curing model, to 
achieve a better prediction of the residual stresses and distortion than the 
approach presented by Svanberg, but with greater computational costs. 
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Warpage optimization was addressed by Gao and Wang [85]. Their approach 
consists of a design of experiments (DOE) to create training points for a meta- 
model, which is then used to link warpage and process parameter. Since they 
focused on injection molding, they used Moldflow during the DOE to create 
the necessary sample points for the meta-model. Thereby, a Kriging model 
was used to link the resulting warpage with the design parameters injection 
time, mold and melt temperature, packing time and pressure. Oliaei et al. [86] 
present a method for the optimization of warpage and volumetric shrinkage, 
utilizing Moldflow and an artificial neural network. Thereby they used a simi- 
lar parameter set as Gao and Wang. 


The curing model presented by Kamal showed that it is able to sufficiently 
predict the resin’s curing behavior, with a reasonable amount of fitting pa- 
rameters and experimental effort. Furthermore, the approach presented by 
Svanberg is an efficient and robust method for the prediction of the occurring 
warpage. As for the draping simulation, computational efficiency is preferred 
over accuracy in this work, since the optimization workflow requires numer- 
ous simulation runs. A detailed process simulation has to be performed after 
the optimization to account for further process influences. 


2.7. Robustness studies in context of optimization 


The idea of incorporating uncertainties from manufacturing in the product 
development process trace back to Taguchi [87, 88], where the topic is cov- 
ered from the final product quality point of view. The first integration of a 
robustness criteria in an optimization workflow was proposed by Conceição 
[89], where they used a robustness measurement as a constraint for a 
weighted sum optimization. 


The first implementation of a robustness measurement as an objective in sin- 
gle-objective optimization was presented by Jin and Sendhoff [90]. They pro- 
posed an expectation-based and a variance-based method for the prediction 
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of the robustness. For the calculation of the robustness value with the expec- 
tation-based method, a number of random sample points in the neighbor- 
hood of the current solution is created, and their mean value is used as a fit- 
ness metric of the current design. The variance-based method utilizes a 
standard deviation around the current solution as an estimation of its robust- 
ness. The standard deviation is calculated a-priori for different conditions, 
based on a reference set. 


For a multi-objective optimization, the first implementation of a robustness 
consideration was presented by Deb and Gupta [91]. In their work, they pre- 
sented two different approaches. The first method is to consider an effective 
fitness function, rather than the real fitness function. Therefore, they pro- 
posed to include deviations in manufacturing or other sources into the fitness 
evaluation step. The second method proposed by Deb and Gupta, is to intro- 
duce a robustness-based fitness function, which is then used as a constraint 
for the optimization. Sun et al. [92] included a variation of relevant design 
parameters in a crash simulation. Utilizing a particle swarm optimization 
method (PSO), they included surrogate models in the optimization loop, to 
avoid computationally expensive simulation runs. 


Several benchmark problems for robustness in optimization runs are pro- 
posed by Gaspar-Cunha et al. [93]. In their work, they included the robustness 
as an objective for the optimization. The robustness of a solution is thereby 
evaluated by a number of test points, randomly created within the neighbor- 
hood of the current solution. An overview of applicable robustness optimiza- 
tion methods was given by Jin and Branke [94]. Furthermore, Beyer and Sen- 
dhoff [95] presented an extensive overview on robust optimization, in which 
they focus on robustness in terms of changes in the production process, 
model sensitivities, and parameter drifts during operation time. A brief survey 
on the consideration of uncertainties in optimization was given by Schuéller 
and Jensen [96], where they suggest to include an uncertainty analysis in the 
product development process to improve the performance and quality of a 
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developed product. A method for the consideration of manufacturing inaccu- 
racies in the topology optimization, to improve the robustness of a solution, 
has been proposed by Troll [97]. Thereby, the focus was on the SMC process 
with varying fiber orientations in the SMC material, but without consideration 
of local reinforcements. 


Since it is important for engineering applications to consider the robustness 
of an optimization result, the robustness should be included as an additional 
decision making criteria for the design engineer. In the context of this work, 
the robustness evaluation is performed during a post processing step of the 
optimization to further account for manufacturing influences. Thereby, the 
aim of the evaluation is to predict a process range, which has to be main- 
tained to ensure that the optimum predicted by the optimization is an opti- 
mum that can be achieved by the manufacturing process. 


2.8. Meta-models in the context of optimization 


In general, a meta-model is used as an approximation of a detailed physical- 
based simulation. It is more efficient to run, which means, that meta-models 
have a trade-off between accuracy and efficiency. Therefore, in context of an 
optimization, a meta-model builds the link between the design parameter 
vector u and the fitness F of the design. If the simulation model s(u) is de- 
fined as 


F = s(u) (6) 


then the meta-model m(u) can be described as 


F = m(u) (7) 
and 


F=Fte, (8) 


where € represents the error of approximation [98]. 
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According to Simpson et al. [98], meta-models are especially useful in early 
states of product design, since there is large uncertainty regarding the con- 
straints, design parameters, and manufacturing process. In addition, meta- 
models are useful for optimization-based robustness studies, since produc- 
tion failures cannot be predicted precisely in an early design state and have 
to be estimated, based on expert knowledge. Simpson et al. [98] give a good 
overview of available meta modeling strategies. 


An extremely flexible but complex modeling approach is the one proposed by 
Krige [98, 99], which is especially well-suited for deterministic problems. 
Thereby, the approach from Krige consists of a known approximation function 
to model the global behavior, and a realization of a stochastic process to cor- 
rect the trend function and ensure the interpolation of all design points. The 
model proposed by Krige is referred to as Kriging model, and was intentionally 
used to create maps of underground geologic deposits using irregular spaced 
boreholes data. A meta-model based optimization approach, for engineering 
applications was presented by Park and Dang [100], where they used an FE 
simulation to create training points for a meta-model, which is then used for 
the optimization. The approach was extended by Rouhi et al. [101], where 
they implemented a loop, with several meta-model creation steps. Thereby, 
the bounds of the sampling space for the meta-model are tighten to the cur- 
rent optimum, to achieve a better representation of the most promising area. 
This is achieved by an refinement of constraints for the sample point selec- 
tion. Furthermore, a good overview on meta modeling in an engineering con- 
text was given by Wang and Shan [102], where they focused on the topics of 
real world applications, role of meta-models in engineering, and applied 
model types. 


Due to its flexibility, the Kriging model is used as meta-model for the robust- 
ness evaluation within this thesis. Furthermore, this type of model provides a 
high prediction quality in areas close to the training points. This characteristic 
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is beneficial in this work, since the robustness evaluation focuses on the vi- 
cinity areas. Thereby, the results gained by the evolutionary algorithm during 
the optimization can be used as training data. 
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3.1. Definition of patch optimization problem 
Some of the content in this section has been published in [103]. 


A general definition of a multi objective optimization problem is given by 
equation (1) in Section 2.4. Based on this the patch optimization problem can 
be formulated as 


Minimize F(u) = [F,(u), F,(w)]", u € U, (9) 


where F; and F, are the objective functions, U represents the design space 
and u the vector of design variables, given by 


u= [x,y,L@]". (10) 


Here, the position of a patch center point M is defined by the set of coordi- 
nates (x,y). Furthermore, @ represents the orientation of the patch, and l 
the patch length. In the proposed approach the width of the patch w is kept 
constant according to usual manufacturing conditions, and is therefore not 
an optimization parameter. The described patch definition is visualized in Fig- 
ure 3-1. 
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Initial component 


Figure 3-1: Visualization of the patch optimization problem, with the center point M, the cen- 
ter point coordinates (x, y), patch orientation , length L, and width w 


Furthermore, constraints for the design variables have to be defined. Based 
on equation (2) and (3), the patch optimization problem is constrained by 


Amin < Xi < Xmax 


Ymin < Yi < Ymax 
(11) 
Emin < li < Imax 


Pmin < Pi < Pmax - 


Thereby the design variables l and @ are modelled as discrete design param- 
eters, with the step size lstep and Pgtep, which restricts the design space to 
practical solutions. 


3.2. Basic definition on the optimization space 


A definition of the optimization space is given to clarify the difference be- 
tween the search, design, and solution space. These three domains are used 
for the definition of the optimization workflow, as well as for the definition of 
the robustness. The basic correlations are visualized in Figure 3-2. Here, the 
design space D corresponds to a set of all possible designs, represented by all 
possible combinations of design variables u, and the solution space S to a set 
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of all possible solutions, represented by all possible fitness values. For an un- 
constrained problem, is 


D,SER. (12) 


The subset F corresponds to the fitness values of all solutions, that are feasi- 
ble for the optimization problem, constrained by the specifications from 
equation (11). Furthermore, the set of feasible designs can be defined by 


U = {ue D|f (u) E€ F}, (13) 


where f(u) is the function used to map between design and solution space, 
which is in the context of this thesis equal to the fitness evaluation. 


Search space 


(o) 


Design space Solution space 


i 


k fw F 


Figure 3-2: Definition of search space, design space, and solution space for an optimization 
problem 


In addition, a search space G is defined. The search space is the representa- 
tion of the design values, used by the optimization algorithm to perform the 
search operations. Therefore a mapping between the search and the design 
space is necessary, where g(u) is a function used to map the design param- 
eter, both ways, between design and search space. In the presented approach 
a binary representation is used for the search space, cf. Chapter 3.3, therefore 
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g(u) represents a transformation from integer values to a binary represen- 
tation. Furthermore, it should be clarified that for optimization approaches 
with real-value-search methods this would lead to the simplification G = D. 


3.3. Evolutionary algorithm 
Some of the content in this section was already published in [103]. 


The idea of an optimization algorithm based on Darwin’s principal of the sur- 
vival of the fittest, was first introduced by Holland [104]. Therefrom, the field 
of evolutionary algorithms evolved, which consists of several different strat- 
egies, like genetic algorithms, evolutionary strategies, genetic programming, 
and evolutionary programming [105]. Today, the most popular type of evolu- 
tionary algorithm is the genetic algorithm, which is why both terms are often 
used synonymously. Contrary to other evolutionary algorithms (e.g. evolu- 
tionary strategies), for genetic algorithms mutation is not the primary search 
operator, but a crossover step is included to propose solutions based on prior 
solutions [106]. Since a genetic algorithm is applied here, this section will only 
focus on this methods. 


The general workflow of a genetic algorithm is shown in Figure 3-3. First, an 
initial population Winitia consisting of W individuals, is created randomly 
within the design space. In the context of evolutionary algorithms, an individ- 
ual corresponds to a design vector u. After a fitness value is assigned for each 
individual (cf. Section 3.4), a selection (cf. Section 3.3.3) and a crossover step 
(cf. Section 3.3.1) is performed to create an offspring generation. In addition 
to the crossover operator, a mutation step (cf. Section 3.3.2) is applied, to 
increase the diversity in the population. 


Once the fitness calculation for the offsprings is done, the new population 
Wi41 is formed by a selection procedure. Thereby, a distinction is made be- 
tween the generational and the steady-state approach. The generational ap- 
proach considers only the offsprings k; for the next generation, while the 
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steady-state approach considers both, parents yw; and offsprings x; [107]. 
Here, the steady-state approach is applied, since it is more likely to converge. 


After the new population w;,, is formed, again a crossover is performed to 
create new individuals. This procedure is repeated until either a defined num- 
ber of generations n is reached, or a convergence criteria is met, cf. Section 
6.3. The result of the optimization is given by the final population Wena), which 
represents the Pareto optimal front. 


In addition it should be mentioned that the size of the population has an im- 
portant influence on the optimizations convergence behavior, since a too 
small population has not enough diversity for a crossover to be useful, while 
a too large population will have a slower convergence rate. 


Initial 
population Offspring Offspring 
(initial) (ki) (ki) 


wy Fitness evaluation „ Mutation . K 
and crossover 
| | 


Fitness 


evaluation 
Crossover | 
ÌE Selection Selection 
«4 7, =: = final 
New Intermediate Final 
population population population 
(isa) (pi + Kj) (rina) 


Figure 3-3: General workflow of an evolutionary algorithm with the main steps fitness evalua- 
tion, crossover, mutation, and selection. 


There is a large variety of different algorithms proposed in literature [48, 108- 
111]. Since the focus of this thesis is not the development of an optimization 
algorithm, but an optimization strategy applying existing algorithms, one al- 
gorithm is selected and used to perform all calculations. The algorithm used 


37 


3. Patch optimization method 


within this thesis is the no-dominated sorting genetic algorithm II (NSGA-II) 
developed by Deb et al. [112], which is an improved version of the NSGA, pro- 
posed by Deb and Agrawal [113]. The peculiarity of the NSGA-II is that the 
fitness calculation is based on a domination value, cf. Section 3.4. 


3.3.1. Crossover 


The purpose of the crossover phase is to create new possible solutions (off- 
spring) based on existing individuals (parents). In general it can be distin- 
guished between real-value and binary crossover operator [113]. For real- 
value operators, the search space is equal to the design space, cf. Figure 3-2. 
In contrast, binary operators use a mapping function g(u) to transfer the 
real-value design parameters into binary representations and, thus, use a bi- 
nary search space for the crossover, cf. Figure 3-2. A binary representation is 
used in the proposed algorithm, since it is the most common approach in lit- 
erature [114]. Therefore, the crossover operators used within the scope of 
this thesis are from the binary type, and only this type will be further dis- 
cussed. 


The first crossover type, which is investigated in this study, is the multi-point 
parameterized binary, where each design variable is threaten individually. 
Therefore a binary representation of the design vector u is created, and a bit 
switching operation is performed on each parameter individually, cf. Figure 
3-4. Here, two parents are used to create two children. Thereby, the parents 
are selected randomly within the current population, and the crossover pro- 
cess is repeated until the desired population size is reached. The behavior of 
the crossover operator can be controlled with the two parameters, probabil- 
ity 9. and number of points u. The parameter u is used to set the number of 
bit-switching operations done per parameter during the crossover process. In 
Figure 3-4 an example with u = 1 is given. The probability parameter Ve is 
used to set the crossover probability for each parameter, and is defined by 
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Ve E [0,1], (14) 


where O represents a crossover probability of 0 %, while 1 corresponds to 
100 %. The crossover probability is used to ensure, that an offspring can also 
consist of unchanged parts of the parents, which leads to a more pronounced 
local search around the parents. 


Parents Offspring 
{x Vi Qi (Xi+1 Vier Pirı lisa} 


i} 
Uia t i L on {0011 001 1000 1010} wisi4 
uiz {1011 M100 1000 101 {1001 M100 100 1010} W412 
Figure 3-4: Example for a multi-point parametrized binary crossover with number of 
points u = 1 


The second crossover operator is the multi-point binary. Here, all design pa- 
rameters are combined in one gene representation, cf. Figure 3-5. The bit- 
switching operation is then performed on this gene, which has the peculiarity 
that one design parameter can be changed on multiple positions, while other 
design parameters may not changed at all. The crossover process is controlled 
with the same parameters, u and Ve, as the multi-point parametrized crosso- 


ver. 

Parents Offspring 

{x vi Gi li} {isa Vier Pirı li+1) 
Uia an 0} {0101110111001010} Ui+1,1 
u;, {1101110010001011} {1101100010001011} u;412 


Figure 3-5: Example for a multi-point binary crossover with number of points u = 2 
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Both binary crossover methods share the characteristic, that the average 
value of the decoded gene stays the same during the crossover process, since 
the bits are purely exchanged cf. Figure 3-6, [113]. This can be proven by tak- 
ing the mapping equation 


m 
giu) = > bir * 2“, 
k=0 


where i = parent1, parent2, offspring1, offspring2 


(15) 
where m represents the maximum number of bits used for the representation 
and the bit-string b; is defined as 


b; = [bo, by, bg], with b, E {0,1}. (16) 


Jparent1 + Iparent Zi Joffspring a+ Joffspring 2 (17) 
2 u 2 i 


The mean value of parents and offsprings is defined by 


With equation (15), and the fact that the number of “ones” and “zeros” has 
to be equal before and after the crossover, it can be seen that the mean value 
has to stay constant during the crossover process. 


Parents Decoded Offspring Decoded 

0110101 61 0110101 119 

010001 33 => goimo0: 25 
Avg. 72 Avg. 72 


Figure 3-6: A special behavior of the binary crossover operator is that the decoded average of 
both parent and offspring strings stays the same 
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In addition, the definition of the crossover spread factor s, taken from Deb 
and Agrawal [113], is introduced by 


= Iparenti — Iparent 2 (18) 


J 
Joffspring 1 — Joffspring 2 


with 
s < 1 are contracting crossover (local search), 


s > 1 are expanding crossover (global search), (19) 


s = 1 stationary crossover. 


Since the mean value is equal before and after crossover, cf. equation (17), a 
stationary crossover means that parents and offspring remain the same be- 
fore and after crossover, which of course is undesirable. 


3.3.2. Mutation 


Beside the crossover, a mutation step is performed as a secondary search op- 
erator. While the primary driver of the optimization is the crossover, the mu- 
tation is used to increase the diversity within the population, cf. Figure 2-7. 
In general, mutation is applied after the crossover step on the offspring pop- 
ulation k;, cf. Figure 3-3. Similar to the crossover, a probability for the muta- 
tion can be defined as 


Om E [0,1]. (20) 


The mutation operation is applied on the binary representation, and switches 
one or more randomly selected bits from O to 1 or vice versa. The process is 
illustrated by example in Figure 3-7. 
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{xi u a  } 
u; {0101 001 1100 1010} 


Ui mutated {0101 0001 1100 1010} 


Figure 3-7: Example for a bit switching mutation operation 


The mutation probability Vm is applied for each parameter individually. There- 
fore, by using basic equations from the theory of probability, the overall prob- 
ability distribution for mutation P can be calculated by 


Pl, Nnm) = (927 * (1-9,)mm) * k (21) 
with 
+ ! 
T- (Nm + Nnm) (22) 
Nnm! * Nm! 


where Nm is the number of mutated and Nam the number of non-mutated 
design variables. 


3.3.3. Selection 


The selection phase is used to choose the individuals that should survive in 
the next generation. Several different strategies are presented in literature, 
like tournament selection, fitness proportional selection, generational selec- 
tion, roulette-wheel selection, elitist selection, and below-limit selection 
[107, 115, 116]. As for the crossover operator, the most common approaches 
are used for the selection. Therefore, the elitist and below-limit selection 
methods are used in this thesis, and this chapter focuses only on these two 
approaches. 
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When using the elitist method, only the individuals on the Pareto front survive 
in the next generation. This means that in the example, given in Figure 3-8, 
only individuals from front 1 are taken from one generation to the next. Con- 
sidering the fronts from Figure 3-8 (a), the elitist approach would not consider 
individual q, when setting up the next generation, since fitness 1 of individual 
qı is slightly better than those of q2. While doing so, promising individuals 
might be rejected when setting up the next generation. To overcome this 
problem, the below-limit approach has been proposed, where an adjusted 
number of fronts is taken from generation to generation. 


While the below-limit approach would work better for the configuration pre- 
sented in Figure 3-8 (a), the contrary applies for Figure 3-8 (b). Here, also the 
individuals on front 2 are preserved to the next generation. Since, all individ- 
uals in a generation serve as possible parents during the crossover, cf. Section 
3.3.1, the convergence rate would be expected to be significantly lower, as- 
suming that parents with considerably worse fitness will result in offsprings 
with worse fitness values. 


(a) 10- (b) 10- 
front 2 front 3 front 2 


064 front 1 front 1 


Fitness 2 


Fitness 2 
o 
> 


0.24 0.24 


T T T T 1 T T T T 1 
0.0 0.2 0.4 06 08 10 0.0 0.2 04 06 08 1.0 
Fitness 1 Fitness 1 


Figure 3-8: Subdivision of the results into different fronts, with fronts good (a) and bad (b) for 
a below-limit selection method 


Both methods can be used in combination with a niching pressure. Thereby, 
only one individual in a defined area is taken. The aim of this method is to 
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increase the diversity along the final front and to prevent local accumulations 
of individuals, cf. Figure 2-7. 


3.4. Fitness calculation 
Some of the content in this section was already published in [103]. 


In terms of optimization, a fitness of a solution is defined as a measure for the 
solution’s capability to fulfill the optimizations objectives [110]. The most 
common approach to combine multiple objectives is to use a weighting func- 
tion, cf. equation (5). However, the weighting approach has the drawback, 
that the optimizations result is directly influenced by the selected weighting 
factors. Therefore, an approach without weighting is used here. 


First, the fitness values of the objective functions have to be evaluated for 
each individual j in the current generation, cf. Figure 3-3. In this work, the 
fitness evaluation is done by an Abaqus simulation, cf. Section 6.2. Once the 
objective functions are evaluated, they must be made comparable. There- 
fore, a domination value dj is calculated 


n 
i=1 


with 


d; 


= f ‚Fı(u;) > F,(u;) and F,(uj) > Fz (uj) 
ji = 


24 
0,else 24) 


where the fitness values F«(u;) represent minimization objectives, n is the 


total number of individuals and 


j = [1,2,...,n]. (25) 
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The selection process is then made, based on the domination value. Figure 
3-9 shows the evaluated domination values for an example set of individuals. 
In addition, the dominated area for each point is highlighted, which corre- 
sponds to the visualization of equation (23). 


1.05 
l 
; 0 ~ 
5 
0.8- . 
2 
a 0.6-4 . 
2 0 
® 4 I 
£ 2 
L 044 
] . 
0.24 0 
Booo 
] 0 
0.04 T T T T 1 
0.0 0.2 0.4 0.6 0.8 1.0 
Fitness 1 


Figure 3-9: Visualization of the domination values and dominated areas for an example set of 
individuals 


Furthermore, the domination value is used to arrange the individuals into dif- 
ferent fronts, cf. Figure 3-8. It is obvious, that all individuals with a domination 
value of zero are those which lie on the Pareto front [46, 47]. 
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4.1. Kinematic draping method in the context of an 
optimization strategy 


Most of the content in this chapter has been published in [103]. 


Due to the high number of simulations performed during an evolutionary op- 
timization, a very efficient draping algorithm is needed. Since kinematic drap- 
ing simulation requires small computation times than constitutive-based 
draping simulation, kinematic methods are, from a computational point of 
view, more suitable for optimization workflows. Furthermore, the optimiza- 
tion shall work as a suggestion for the designer and not as a final verification. 
Hence, a mechanical draping simulation must be performed afterwards to en- 
sure manufacturability of the optimized patch setup without forming defects 
like wrinkling, fiber fracture or gapping. This is similar to topology optimiza- 
tion, where a verification simulation is needed to ensure the load-bearing ca- 
pacity of the optimized topology [7]. The focus of the draping simulation 
within the presented optimization workflow is the calculation of the final 
patch geometry, while the prediction of forming defects is not yet relevant. 
Therefore, kinematic draping simulation is chosen to be sufficient for the as- 
pired purpose. The drawbacks regarding the dependence on the initial start- 
ing point are not highly important, since all possible starting points are taken 
into account by the proposed optimization workflow. Furthermore, a com- 
parison of FE based and kinematic draping simulation will be given in Section 
4.3 using the example of a curved reference structure. This comparison is 
used to discuss the advantages and disadvantages of both methods. 
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4.2. Workflow of the draping simulation 
The method presented in this chapter was already published in [103]. 


The general assumption of kinematic draping simulation is the inextensibility 
of the material in fiber direction. Depending on the type of fiber reinforce- 
ment, the deformation of the material is assumed as pure shear or simple 
shear along the fiber direction. While woven fabrics deform in pure shear, 
unidirectional reinforcements are assumed to deform in simple shear, as 
shown in Figure 2-8. 


Starting point a, patch orientation @, patch length lọ and width wọ are the 
necessary input parameters for the kinematic draping simulation. Another re- 
quired input is the geometry of the initial component, on which the patch is 
going to be draped. The first step is to create the initial paths A and B on the 
part surface, starting from the given starting point in the center of the patch, 
cf. Figure 4-1 (a). Thereby, path A is the initial path in fiber direction, while B 
is perpendicular to the fiber direction. 

For the calculation of the next node a<'*1> from the current node a<'> at 
path A, a number of constraints need to be fulfilled. First, the distance be- 
tween both nodes is set to be equal to the defined mesh size m 


<i+1> _ „<i> (26) 


a a =m. 


Furthermore, the node a<'*1> has to be on the initial component surface 
ate = Fas’, as), (27) 


where axt?>, a$!*!>and a$'*1*are the x, y, and z components of the new 
node and F(x, y) represents the surface of the initial component. Addition- 
ally, it has to be ensured that the given global patch direction @ is maintained. 
With these assumptions all nodes along the initial paths A and B can be cal- 


culated, cf. Figure 4-1 (a). 
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asAOBU+1)> 


Figure 4-1: Node calculation, starting from the initial paths A and B (a), and for all further 
nodes (b). 


Based on these initial paths, all other nodes can be calculated. From the as- 
sumption of simple shear (Figure 2-8) the following equations can be derived 
for the parallel paths 


jase neg u aA Da(>| =m, (28) 


ERDE? _ aa >| =m, (29) 


Io wo 
with i = 1,2,...—, J =1,2,..— 
m m 


Equation (28) and (29) represent that the distance between the fibers as well 


as the fiber length is constant during the draping process. Here, a<4+).8U)> 


<A(,BU+)> is the closest node on path B, 


is the next node on path A, while a 
cf. Figure 4-1 (b), lọ the initial patch length and wọ the initial patch width. 
Again, the node has to be on the part surface, which is represented by Equa- 
tion (27). Rearrangement of equations (27), (28) and (29) leads to the follow- 


ing system of equations 
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which has to be solved for all nodes. The presented kinematic draping ap- 
proach is implemented in MatLab and embedded in the optimization work- 
flow. 


4.3. Numerical verification 


To demonstrate the ability of the kinematic draping method used within this 
thesis, it is compared with results from an established mechanical method 
developed by Dorr et al. [117]. Therefore, three reference patch configura- 
tions are defined, cf. Figure 4-2. The tool geometry used here has different 
angles for the two edges D1 and D2. 


Configuration A is a simple, single-curved forming of the patch at the edges. 
In configuration B the patch is placed over the corner bendings, which results 
in more complex deformation than configuration A. In the last configuration, 
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the patch is additionally rotated to further increase the complexity of the 
forming process. 


A [0/90/0] layup with a total thickness of 1.5 mm is used for the mechanical 
draping simulation. For each reference configuration, the first 0° layer is com- 
pared with the results from the kinematic draping simulation. 


Figure 4-2: Tool geometry with two different angles at the sides D1 and D2, and three refer- 
ence patch configurations A, B, C 


The results of configuration A are shown in Figure 4-3. For the simple defor- 
mation applied here, kinematic and mechanical draping simulation create the 
same results. This shows, that the kinematic draping is able to create reliable 
results for simple geometries. 
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Figure 4-3: Comparison of the results for configuration A: kinematic method (green) is congru- 
ent with mechanical method (blue) [118]) 


Comparing the results from configuration B, cf. Figure 4-4, it can be seen that 
the kinematic draping method is not able to give a correct representation of 
the deformation process within the corner bendings. The corner bendings re- 
sult in material accumulations, which cannot be modelled by kinematic drap- 
ing methods, cf. Figure 4-4 (b). Therefore, the folding of the patch in the cor- 
ner bending can not be predicted correctly. The disadvantages of the 
kinematic draping simulation, when it comes to doubled curved structures, 
has also been demonstrated by Sharma and Sutcliffe [70]. Comparing the re- 
sults of the simulations, cf. Figure 4-4 (c), with the experimental result, cf. 
Figure 4-4 (d), it can be seen that the mechanical method predicts the final 
geometry better than the kinematic draping simulation. However, the main 
part of the patch was draped correctly, while the differences are only in small 
areas of the patch. Therefore, the result of the kinematic draping simulation 
can be seen as a good approximation of the patch’s final position. 
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Figure 4-4: Comparison of the results for configuration B: the kinematic model (green) differs 
from the mechanical model (blue) only in the corners of the structure, perspective 
view (a), top (b), bottom view (c) [118], experimental result (d) [119] 


The results of the last configuration are presented in Figure 4-5. As for config- 
uration B, the kinematic draping simulation shows disadvantages in predict- 
ing the correct geometry in the corner bendings. Here, the kinematic draping 
simulation is not able to correctly predict the material movement resulting 
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from the corner bending. However, as for configuration B, the differences af- 
fect only small areas of the patch, while the main parts of the geometry are 
predicted correctly. Furthermore, it should be mentioned that this effect de- 
pends strongly on the considered geometry, since the deviation created at 
the double curved geometry would increase with increasing patch length be- 
yond the curvature 


Figure 4-5: Comparison of the results for configuration C: top view with details D1 and D2 
(green: kinematic method blue: mechanical method [118]) 
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4.4. Conclusion on draping simulation 


The advantage of the kinematic draping simulation becomes clear, when 
comparing the necessary time for the simulation, cf. Table 4-1. While the kin- 
ematic draping simulation is performed within only a few seconds, the me- 
chanical method needs a couple of hours. Since a large number of draping 
simulations is necessary during one optimization run, it becomes obvious that 
a mechanical method cannot be used here. 


Table 4-1: Comparison of the computational times for the reference configurations 


Configuration Kinematic draping Mechanical method 
A 25s 1h54 min 
B 25s 1h 23 min 
C 30s 8h 31 min 


4.4. Conclusion on draping simulation 


A kinematic draping method has been implemented, which fulfills the re- 
quirements of an optimization strategy based on an evolutionary algorithm. 
A comparison of kinematic and mechanical draping methods has been pre- 
sented to show the advantages and disadvantages of a kinematic draping ap- 
proach. Results and computational time have been compared for three dif- 
ferent reference patch configurations, with different degrees of geometrical 
complexity. Taking the results of the comparison with the mechanical draping 
simulation into account, the kinematic draping simulation can be seen as a 
good approximation of the patch’s forming behavior. In general, for the kine- 
matic draping simulation, no material parameters are necessary. In addition, 
the computational time of a kinematic draping simulation is much faster, 
compared with the mechanical method. 


Furthermore, the optimization shall work as a suggestion for the designer and 
not as a final verification. Therefore, a mechanical draping simulation must 
be performed afterwards to ensure manufacturability of the optimized patch 
setup without forming defects like wrinkling, fiber fracture or gapping. 
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5 Curing and warpage simulation 


5.1. Kamal-Malkin curing model 


The calculation of the warpage, occurring during the manufacturing process 
is directly related to the curing of the resin. Therefore, the approach used for 
the prediction of the curing is presented first. For an exothermal crosslinking 
reaction, the degree of cure is defined by 


H(t) 
H max 


q(t) = (31) 
where H(t) and Hmax are the accumulated enthalpy until time t and the total 
reaction enthalpy [84]. In addition, q(t) is per definition 


qe [0,1], (32) 


where 0 represents the uncured and 1 the fully cured state. Based on this, a 
cure rate 


d 
I fT), (33) 


dt 
can be defined depending on the state of cure q and the temperature T, 
where the function f (q, T) depends on the resin type. Several different ap- 
proaches to describe the cure rate, depending on the used material type, can 
be found in literature. An overview is given by Halley and Mackay [76]. 


The simplest form to model the development of cure is the n order 
model [76] 
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== ka-q)" (34) 


with the Arrhenius equation 


k = kọ exp (- =) Ba) 


with R as the universal gas constant, the material dependent parameter ky 
and A, and the order of reaction n. In addition, an autocatalytic model ap- 
proach can be given with the kinetic equation [74, 120] 


t = u 36 
dt 2. kq"(1 q)” , ( ) 


where k is defined according to equation (35), while m and n represent the 
order of reaction. Based on this, Kamal [73] combined both approaches and 
derived the following phenomenological model 


dq (37) 
dt = (kı +k,q™)(1 - q)” . 
Again, k; is calculated with the Arrhenius equation 
AN: 
ky = Kig €XD (- =) TEIR (38) 


Since this model is used in this work to describe the curing behavior, the six 
material dependent parameters n,m, kio, k20, A1, and A, must be deter- 
mined. The data have been determined by DSC measurements performed by 
Schwab [121]. 


It should be mentioned that all of these approaches have the drawback, that 
a small initial state of cure has to be defined to start the reaction. Therefore, 
the initial state of cure is set to 


qo = 0.01. (39) 
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5.2. Svanberg method for expansional strain 


5.2. Svanberg method for expansional strain 


The expansional strain is calculated based on the approach proposed by 
Svanberg and Holmberg [80, 82]. In addition, for the calculation of cure-in- 
duced residual stresses and deformation in fiber-reinforced composites, ei- 
ther the material properties have to be homogenized or the fiber and matrix 
have to be considered separately. Hence, homogenized material parameter 
will be used here, the mapping and homogenization step is explained in Sec- 
tion 5.4. 


The stress tensor is defined by 
Gij = ijet Ekt — Cijkı Eka (40) 


where Cjjx;, is the stiffness tensor, €x, is the mechanical and ef the expan- 
sional strain tensor. Thereby, the stiffness tensor is temperature dependent, 


no stress ,q < ge, and T = T,(q) 
Ciji = $ Cijev 2 Ageı and T > T,(q) (41) 
Ciri ‚T < T,(q) 


where Cir and Chr are the stiffness tensor for the rubbery and glassy state, 
while in the liquid state no stress development is expected to occur. Consid- 
ering the approach from Svanberg and Holmberg [80, 82], the expansional 
strain can be divided into a thermal part £f, and a chemical part e£, 


Ein = Ea + Ex (42) 
The thermal strains can be expressed by 
(43) 


T oT 1 
Ek = | aard ae , 


0 
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5. Curing and warpage simulation 


where the coefficient of thermal expansion a,; depends on temperature and 
degree of cure 


ak ‚q < Agel and T > T,(q) 
ari = \ Ca q 2 ge and T > T,(q) (44) 


where al,, at,, and akı are the coefficients of thermal expansion in the liquid, 
rubbery, and glassy state. T, is the glass transition temperature and qge the 
gel point of the resin. The gel point characterizes the change from liquid to 
solid behavior. 


Similar to the coefficient of thermal expansion, a coefficient of chemical 
shrinkage is introduced, which correlates the chemical expansion to the cur- 
ing rate. Thus, the strain resulting from the chemical shrinkage can be written 
as 


t 
ðq, 
& = [Buero hat , (45) 
0 


while the coefficient of chemical shrinkage Fk is also dependent on the tem- 
perature and on the degree of cure, and can be expressed by 


Bra ds Agel and T > T,(q) 
Pri = Bra 2 Agel and T > T(q) (46) 
Bii T < Tg(q) 


where ß},, Bi), and BE, are the coefficients of chemical shrinkage in the liquid, 
rubbery, and glassy state. Since the chemical shrinkage is only present in the 
matrix, the coefficient Pı can be calculated from the matrix properties. 
Therefore, Svanberg and Holmberg [80, 82] suggested the following approach 
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5.2. Svanberg method for expansional strain 


i Ay matrix 
Ban ies ; (47) 
3 
where AV™atrix represents the relative matrix cure shrinkage, which repre- 
sents the relative volume change of a resin from liquid to fully cured state. 
This corresponds to the assumption, that the cure shrinkage is evenly distrib- 
uted in all spatial directions. Furthermore, the following simplification shall 
apply 
matrix _ pmatrix _ pmatrix 
XX = byy — Pzz 
matrix =0 (48) 
with k,l = [x,y,z],k #1 


During the homogenization step, the coefficients yx, Pyy, and pzz are ad- 
justed according to the corresponding fiber volume fraction and fiber orien- 
tation, cf. Section 5.4. 


The relative cure shrinkage of the matrix can be estimated from an experi- 
ment by [80, 81] 


liquid d 
matrix _ Pmatrix(Teure) — Pmatrix (Teure) 1 (49) 
Te 
Pmatrix (Teure) q (Teure) 
h liquid cured 5 RR . 
where pnatrix and Pmatrix are the density of the liquid and the cured resin, 


respectively. Furthermore, Teure is the temperature at which the resin was 
cured during the experiment and q(T,„re) is the degree of cure achieved at 
the end of the experiment. 


The glass transition temperature is the temperature at which the transfor- 
mation from rubbery to brittle state occurs, resulting from temperature de- 
crease. Hence, the knowledge of this value is essential for the definition of 
the material properties stiffness, thermal expansion and chemical shrinkage. 
The glass transition takes usually place within a temperature range, rather 
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than at one specific temperature. The glass transition temperature depends 
on the degree of cure, which can be expressed with the DiBenedetto equation 
[80, 82] 


Te —T, x 
eS (50) 
Toe Ten = (l= A) eq 


where Tg, and T,.. are the glass transition temperature of the uncured and 
the fully cured matrix, respectively, and A represents a material system de- 
pendent parameter. 


5.3. Implementation of the curing and warpage 
model 


Svanberg’s approach is implemented within an Abaqus user subroutine of the 
type UEXPAN. To be implemented, equations (42) - (45) have to be formu- 
lated in an incremental form, using the forward integration scheme 


ef (t + At) = cE (t) + Ach, (51) 


where e® represents the expansional strain. According to the approach, the 
expansional strain can be split up in a thermal (£T) and a chemical (£f) com- 


ponent 
Ack, = Ae}, + Act. (52) 
with 
Ae}, = ayı * AT (53) 
and 
Aegi = Bri * Aq, (54) 
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5.4. Calculation of the local material parameter 


where AT and Aq are the temperature and cure gradient. The coefficients a; 
and pk; are selected according to equations (44) and (46). For the used mate- 
rial, the coefficients al, and £}, of the liquid state are assumed to be zero, 
since the strains remain undeveloped in the liquid state. In addition, the cur- 
ing coefficients in rubbery state By, and glassy state Be, are considered to be 
identical. Besides this, the stiffness tensor depends also on temperature and 
degree of cure, cf. equation (41). However, since there is no material data 
available for the used material, the stiffness tensor for the glassy state is used 
for all states. This simplification should be sufficient enough here, since the 
main focus of this thesis is on the optimization procedure. 


Since an initial degree of cure gu is required for the implemented curing 
model, cf. equation (37) and (38), it is set to 0.1 %. Furthermore, the gel point 
Age) is set to 60 % for all calculations performed. Both values represent as- 
sumptions based on known material systems. Since no measurement for the 
development of the glass transition temperature T, dependent on the degree 
of cure is available for the material system used, T, is assumed to be constant 
during the curing process. 


5.4. Calculation of the local material parameter 
based on fiber orientation and fiber volume 
content 


The coefficient of thermal expansion for the patch is calculated based on the 
method proposed by Schapery [122]. According to Schapery’s approach, the 
coefficient of thermal expansion in fiber direction can be calculated by 


a, = matrix * tmatrix * (1 — 9) + Brier * @tiber*@ (gy 
Ematrix = (1 =. 6) + Efiber *0 


where Efiber and Ematix are the Young’s moduls of fiber and matrix, while 8 is 
the fiber volume content. Here, the carbon fiber is assumed to be isotropic. 
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5. Curing and warpage simulation 


This simplification should be feasible, since the focus of this work is on the 
development of an optimization procedure. Furthermore, the coefficient of 
thermal expansion perpendicular to the fiber direction is 


a= (1 + Nmatrix) * Æmatrix * (1 ze 0) + (1 + Nfiber) i (56) 
Qfiber * 0 —- a * (matrix = (1 = 0) + Nfiber * 0), 


where Nmatrix and Nfiber are the Poisson’s ratios of matrix and fiber. 


The fiber orientation distribution within the discontinuous fiber reinforced 
component is calculated by a preceding process simulation. The resulting an- 
isotropic coefficients of thermal expansion are calculated with the method 
proposed by Rosen and Hashin [123]: 


Amatrix — Xfiber 


-1 —1 
Katrix = Kiber 


þess- | 
* * 22 | ———__ 5 
oy Kmatrix Kfiber 


Qij = Amatrix * (1 z 0) + Qfiber * 0+ 
(57) 


where Skrij is the compliance matrix, and Kfiber and Kmatrix are the bulk mod- 
ulus of fiber and matrix. 


The coefficients for chemical shrinkage, also dependent on fiber orientation 
and fiber volume content, are calculated with the same approach as the co- 
efficient of thermal expansion. The material parameters used for the thermal 
and chemical coefficients are given below in Table 5.2. 


5.5. Numerical verification 


Since the curing model is the basis for the calculation of inner stresses and 
warpage, the results obtained with the curing model are investigated first. 
Experimental results are taken from Schwab [121]. For the characterisation 
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5.5. Numerical verification 


of the used matrix system, DSC measurements have been performed by 
Schwab to measure the heat flux created during curing. The total enthalpy 
Hmax and the accumulated enthalpy H(t) of the exothermal crosslinking pro- 
cess can be calculated with a numerical integration. The temperature profiles 
used for the experiments and the simulations are given in Figure 5-1. Thereby, 
two different heating rates are used to characterise the curing kinetic of the 
resin. Based on this, the degree of cure over time can be described by equa- 
tion (31) and the fitting parameters of the Kamal-Malkin model can be calcu- 
lated (equations (37) and (38)), which are given in Table 5-1. 


Table 5-1: Fitting parameters of the Kamal-Malkin curing model, based on DSC measure- 
ments performed by Schwab [121] 


kıo kzo Ay monl A, noll m n 
0.9601 0.9996 3.2637 * 10* 3.0737 * 10° 0.7743 1.3988 
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Figure 5-1: Temperature curve of the heat rates used for DSC experiments and simulations 
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5. Curing and warpage simulation 


The comparison of experimental and simulation results is shown in Figure 5-2. 
Here, cure rate and degree of cure are compared in both the time and tem- 
perature domains. Comparing the results in the time domain, a good agree- 
ment of experiment and simulation can be seen for a heat rate of 20 K/min, 
while for a lower heat rate of 15 K/min the cure rate is slightly overestimated 
by the simulation, cf. Figure 5-2 (a). When looking at the degree of cure, cf. 
Figure 5-2 (b), the first increase is slightly too early for the higher heat rate, 
and a little delayed for the lower heat rate, showing an acceptable overall 
compromise. Comparing the results in the temperature domain, a small shift 
in the cure rates to higher temperatures can be observed for both heat rates, 
cf. Figure 5-2 (c). The degree of cure fits well for the higher cure rate, whereas 
there is a small deviation for the lower cure rate, where the reaction starts at 
a higher temperature in the simulation than in the experiment, cf. Figure 
5-2 (d). In general, only small deviations between experiment and simulation 
occur, which shows that the Kamal-Malkin curing model has a good prediction 
quality for the curing behavior of the material system. 
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Figure 5-2: Comparison of the results from the DSC experiments (dashed lines) and simulations 
(solid lines) for two different heat rates (15 K/min, 20 K/min): development of the 
cure rate (a) and degree of cure (b) over time; development of the cure rate (c) and 
degree of cure (d) over process temperature 


After the curing model, the strain calculation approach is discussed. There- 
fore, two different material model setups are used. The first model considers 
only the different coefficients of thermal expansion for fiber and matrix, while 
the second model uses the complete Svanberg approach, including chemical 
shrinkage, cf. Section 5.2. The comparison of the two models shall demon- 
strate the influence of the chemical curing on the strain and warpage devel- 
opment. The simulation model and the temperature profile used for the eval- 
uation are shown in Figure 5-3, while the material parameters used for the 
calculations are given in Table 5.2. 
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5. Curing and warpage simulation 


For simplicity, a unidirectional fiber direction with a symmetric layup is used, 
with the fiber direction equivalent to the global x-direction, cf. Figure 5-3. The 
temperature profile consists of a period of constant temperature, which sim- 
ulates the time the material spends in the tool, and a cool down phase, in 
which the temperature is reduced to room temperature. The temperature is 
thereby applied uniformly on the whole part. To exclude an influence of geo- 
metrical effects on the results, a flat plate structure is used. To enable a free 
deformation in membrane direction, the nodes B1, B2, and B3 are con- 
strained in x, y, and z direction, respectively. 
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Table 5-2: Material parameters used for the calculation of the expansional strain 


Material parameter Symbol Value Unit Source 
Glass transition tem- Tg 140 °C fitting 
perature 
Initial degree of cur- qo 0.1 % assumption 
ing 
Gel point Agel 60 % assumption 
Volume shrinkage AV 5.6 % [124] 
Young’s modulus E matrix 3.1 GPa [125] 

Efiber 73 GPa [125] 
Poisson’s ratio Nmatrix 0.39 - [125] 

Nfiber 0.22 - [125] 
Fiber volume content smc 0.22 - [125] 

Bitch 0.47 - [125] 
Glassy state 
Thermal expansion aË tix 6 * 10° 1/°C [126] 

afer  225*10 APCs [127] 

6 

Curing P -0.0187 1/°C Equation (47) 

Piber 0 1°C - 
Rubbery state 
Thermal expansion Onatix 1,5*10% IC. 25 Ge. [80] 

Kfiber 2,25*10 1/°C [127] 

6 

Curing Pratis -0.0187 1/°C Equation (47) 


Piber 0 1/°C 7 
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Comparing both simulation results, it can be seen that the developed strain is 
significantly underestimated when curing effects are not considered, cf. Fig- 
ure 5-4 (b). The chemical shrinkage occurs after the gel point is reached (q = 
gei), until the temperature is below T,. Thereafter, no further curing takes 
place and the strain development depends only on the directional-dependent 
coefficient of thermal expansion. Consequently, from there on the strain de- 
velopment for both examples is alike. 


Cooling rate (25K /min)] 


Temperature [°C] 


Time [s} 


Figure 5-3: Plate model with the constraints B1, B2, B3, and the reference evaluation point A 
(a), and temperature profile (b) 


Comparing the strain components, the strain transverse to the fiber direction 
(cf. Figure 5-4 (b), E22) is largely effected by the influence of the chemical 
shrinkage, since the matrix behavior is dominant. Contrary to this, the influ- 
ence is less significant in the fiber direction (cf. Figure 5-4 (b), E11). 
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Figure 5-4: Resulting displacement (a) and resulting strain (b) at reference point A, cf. Figure 
5-3 (a) with consideration of curing and without consideration of curing. 


5.6. Conclusion 


The curing model presented by Kamal [73] was used to model the curing be- 
havior or the UPPH resin system used within this thesis. A comparison of the 
simulation results with experimental results showed a good accordance for 
different heat rates. 


Beside the curing calculation, a strain calculation was performed based on the 
method proposed by Svanberg [82]. With this approach, the influence of the 
incorporation of chemical strain during the curing process could be demon- 
strated. For this purpose, the strain development with consideration of curing 
was compared with a simulation, where only different coefficients of thermal 
expansion where considered. 


With dependence on fiber orientation and fiber volume content, the effective 
coefficients of thermal expansion for the UD reinforcement patch were cal- 
culated with the method of Schapery [122]. The anisotropic coefficients for 
the discontinuous fiber reinforced component were calculated with the ap- 
proach proposed by Rosen and Hashin [123]. Furthermore, the assumption 
was made, that the effective coefficients of chemical shrinkage can be calcu- 
lated with the same approach as the coefficients of thermal expansion. 
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6 Multi-objective optimization along 
the CAE chain 


6.1. Manufacturing constraints 


The incorporation of manufacturing constraints within the optimization gives 
the essential benefit, that the optimized part can be manufactured as pre- 
dicted or at least close to the predicted optimum. Therefore, the manufactur- 
ing constraints, which apply for the SMC process and should be integrated 
within the optimization loop, will be discussed in this section. First, the basic 
constraints of the compression molding process with local reinforcement 
patches (co-molding), and their influence on the optimization are presented. 
Thereafter, variations resulting from the production of the semi-finished 
product and the co-molding process are presented, which could influence the 
optimization result. 


Semi-finished product 


Since the cutting of the patches is done automatically, usually no scattering 
in the initial length and width of the patch is expected [128]. However, the 
positioning of the prepreg material on the cutting table is done manually. As 
a result, the initial fiber orientation may deviate from the patch length direc- 
tion, cf. Figure 6-1 (a). Thereby, the variation of the orientation is typically in 
a range of +3° from the optimum orientation [128]. Additionally, the fiber 
orientation can be influenced by the prepreg production process, in a way 
that a misalignments goes through the whole patch, cf. Figure 6-1 (b). This 
failure can result in rated break point of the patch. The variations caused by 
fiber misorientation are considered by a variation of the orientation of the 
total patch within the robustness evaluation, cf. Section 8.3. 
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Furthermore, to minimize waste and simplify the manufacturing process, the 
reinforcement patch should have an initially rectangular geometry [128]. 
Thereby, a minimum and a maximum patch size, depending on the initial com- 
ponents dimensions, should be considered. The minimum and maximum con- 
straint is directly implemented in the definition of the optimization problem, 
cf. equation (11), and therefore implemented in the optimization algorithm. 


(b) 


Figure 6-1: Inaccuracies during the prepreg production: misalignment of patch and fiber direc- 


tion during cutting (a), and shift in the fiber orientation of the semi-finished 
product (b) 


Draping process 


An obvious restriction is that the reinforcement patch has to be completely 
within the area of the initial component. A further restriction, resulting from 
the process, is that a patch reinforcement is only possible on either one or 
the other side of the initial component. These restrictions are incorporated in 
the draping simulation, where the patches are cut at the border of the initial 
component, cf. Section 4.2. Furthermore, the final geometry of the reinforce- 
ment patch is predicted by the kinematic draping simulation. A further pro- 
cess constraint, incorporated within the draping simulation is the prevention 
of undercuts of the patch, which would be impossible to manufacture. In ad- 
dition, the consideration of user-defined no-go-areas for the reinforcement 
patch is possible within the draping simulation. No-go-areas in the context of 
patch reinforcements, to be considered here, are preserved spaces for inserts 
or areas with too complex geometry or a too high slope, making the draping 
process impossible [128]. 
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Co-molding process 


The third group of influencing variables are the process-induced local varia- 
tions of the material constitution. Therefore, the manufacturing influence re- 
lated to the DiCoFRP component is the fiber orientation, which is provided by 
a process simulation and used as input for the optimization. Thereby, the fiber 
orientation influences the mechanical performance and the warpage, which 
is predicted with respect to the fiber orientation, cf. Section 5.4. 


During the co-molding process no influences regarding the thickness of the 
patch are expected [128]. The deformation and displacement of reinforce- 
ment patches has been studied by Biicheler and Henning [129]. Their results 
for the displacement and deformation of the patch during the molding pro- 
cess are shown in Figure 6-2. Since the material system used within this thesis 
is a UPPH resin, the deviations for the UPPH resin are going to be used as 
input for the robustness evaluation, cf. Section 8.4. Once the co-molding pro- 
cess is done, the adhesion between the patch and the SMC material is un- 
problematic, since the same resin system is used for both materials. 
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Figure 6-2: Resulting total deformation (left) and displacement (right) for a patch reinforce- 
ment during the manufacturing process [129] 


Summarizing the influences of the manufacturing inaccuracies, two main 
sources of deviations can be identified, orientation, in terms of the whole 
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patch or the fibers, and size changes, like patch length and width. Further- 
more, the most important question, to be answered by the robustness eval- 
uation is, how accurate the process has to be, to ensure that a predicted op- 
timum remains an optimum for the real part. Therefore, a robustness 
evaluation is performed after the optimization. 


6.2. Workflow of the multi-objective patch optimi- 
zation 


Most of the content in this chapter was already published in [103]. 


The general approach of the patch optimization chain is shown in Figure 6-3. 
Thereby, the optimization procedure can be divided into the three main parts: 
pre-processing, optimization loop and post-processing. The pre-processing 
phase starts with a design phase for the initial component. After the initial 
design is specified, a process simulation is performed to predict the fiber ori- 
entations for the SMC component. Based on the fiber orientations, the mate- 
rial parameters are calculated within a mapping step. After the pre-processing 
part is done, the optimization starts. Since the optimization loop is the main 
scope of this work, it will be discussed in the following sections. The last part 
of the optimization chain is the post-processing step. Here the Pareto front is 
evaluated (Section 6.3) and heat maps are created, cf. Section 6.4. 
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Figure 6-3: CAE chain used for the patch optimization 
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6.2. Workflow of the multi-objective patch optimization 


The more detailed workflow of the optimization loop is shown in Figure 6-4 
(a). Firstly, the initial population of design variable vectors u; is created with 
uniform random distribution for each parameter within the given design 
space. Secondly, the initial fitness calculation is performed (Figure 6-4 (c)), 
where the kinematic draping simulation (cf. Chapter 1) is integrated into the 
fitness calculation step (cf. Section 3.4) of the multi-objective genetic algo- 
rithm. The draping simulation is used to create the necessary input for the 
structural simulation in terms of position and fiber orientation of the formed 
CoFRP patch. Furthermore, the draping simulation is used to return the fit- 
ness value F,. Fitness value F, represents the amount of used patch length. 
Since the random selection of design variables u may create patch lengths, 
which overlap the boundaries of the geometry, two different scenarios have 
to be distinguished: 


1. The patch fits completely on the DiCoFRP component: Here the 
patch usage is equal to the design parameter, representing the patch 
length. 

2. The patch fits not completely on the DiCoFRP component: Here the 
patch is cut at the border of the component and the resulting length 
is returned as F3. 


Subsequent to draping simulation, a structural simulation computes the fit- 
ness value F,, which represents the compliance of the component, quantified 
by the global strain energy. For linear-elastic material behavior, the strain en- 
ergy can be defined by 


n 
1 
w.ebal = | o'’xedV = =), P} * di (58) 
a i=1 


where P; and d; are the nodal forces and nodal displacements of a finite ele- 
ment model with n element nodes. 
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Figure 6-4: Integration of the fitness calculation in the genetic optimization algorithm of the 


Final 
Population 


Dakota tool box (a), used to solve the patch optimization problem (b) by applying 
the two-step fitness calculation (c). 


The fitness calculation step is slightly adapted, when the warpage simulation 
is included in the optimization run, since the fitness value F, results then from 
the warpage simulation. Therefore, the fitness calculation is performed as 
presented in Figure 6-5. In this case, the draping simulation is only used to 
calculate the final patch geometry and the resulting fiber orientation as input 
for structural and warpage simulation. The resulting patch length is not inves- 
tigated as an additional objective function in order to limit the number of fit- 
ness values to be evaluated. 


Draping Simulation 


Structural Simulation Warpage Simulation 
F; (u) Fz (u) 


Figure 6-5: Fitness calculation when incorporating a warpage simulation in the optimization 
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For the optimization algorithm, the open access software toolkit Dakota is 
used, an implementation of the Sandia National Laboratories [130]. The Da- 
kota toolkit provides diverse operators for genetic algorithms (as the ones 
described in Section 3.3) as well as an interface between analysis code and 
optimization methods. Within the optimization loop, two termination criteria 
can be used: a maximum number of iterations n and a convergence metric. 


6.3. Convergence criteria 


Since multi-objective optimizations have the goal to predict a wide range of 
the Pareto front, more than one criterion is needed to properly characterize 
the performance. The criteria presented first in this section, are developed to 
evaluate the final results of the optimization. In addition, a convergence met- 
ric is necessary during the optimization runs to measure the improvement 
between two generations. Therefore, Dakota offers a metric consisting of 
three different criteria, which cannot be changed manually [130]. Finally the 
criteria developed are compared with the Dakota criteria. 


The three criteria, introduced within this thesis to evaluate the convergence 
rate, the quantity of solutions, and the diversity, are: 


1. Dominated area for the evaluation of the convergence rate. The domi- 
nated area is defined as the area enclosed by the current front and the 
extreme points k, and k, and therefore represents the progress of the 
optimization. Therefore, the extreme points are defined by 


ki = LF, (uj), F2 (u;)|",i € 1,2 (59) 
ky = {k|F; (u1) > F,(u) Vu E€ U} (60) 
kz = {k|F,(uz) > F,(u) Vu € U} (61) 


The dominated areas of two successive fronts are exemplarily shown in 
Figure 6-6. 
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2. Number of solutions along the Pareto front. 

3. Patch length distribution. The patch length is modeled as a discrete de- 
sign variable. Therefore, the percentage share of the patch length parti- 
tions along the Pareto front is used to describe the distribution along the 
front. This criterion is only used, if the patch usage is one of the objec- 
tives. 
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° 
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s 
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N 
uù 
0.24 “ 
| kı 
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0,0 0,2 0,4 0,6 0,8 1,0 


F1 (standardized) 


Figure 6-6: Dominated area calculation for generation i (black solid line) and i + 1 (grey dot- 
ted line), and the maximum points k4 and k, (red dots). 


The following criteria are implemented in the Dakota workflow, and will 
therefore be used during the optimization run: 


1. Spread of the solutions within the design space, which is calculated by 
AF = |k; — k,| (62) 
where kı and kz are calculated according to equation (60) and (61). 
2. The density along the current Pareto front, is used to measure the diver- 


sity within a current solution. 
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3. The relative “goodness” of the current front. To measure the relative 
goodness, the number of individuals on front (i — 1) dominated by front 
i is counted. 


All three criteria are evaluated for each generation and implemented as per- 
centage change between two generations. From these three criteria the high- 
est percentage change is taken to evaluate each generation. The optimization 
is terminated, when a user defined percentage change is less than a specified 
value, for a defined number of subsequent generations. [130] 


Comparing the proposed convergence criteria with the criteria implemented 
in Dakota, a number of advantages can be seen. When evaluating the conver- 
gence rate, the Dominated area is used instead of the solutions spread, since 
it is a more precise method to measure the optimizations progress, taking all 
points on the front into account instead of only the extreme points. The di- 
versity measurement is done with the patch length distribution, since it sub- 
divides the solutions in groups connected to the input parameter. The overall 
progress of the optimization can be better described, when taking the domi- 
nated area and the total number of solutions into account, instead of the Da- 
kota-defined “goodness” of a generation. 


6.4. Heat-maps 


Since the multi-objective optimization does not result in one best solution, 
but in a set of solutions, a tool needs to be introduced to help a designer to 
choose the right position for the patch reinforcement. The final Pareto front 
can be used to select one specific solution, but thereby a large amount of the 
information gained during the optimization, is not used any more. Therefore, 
the heat map is proposed here as an additional tool to visualize the optimiza- 
tion results. 


The heat map focuses on the evaluation of the global strain objective in con- 
junction with the patch usage or the warpage objective. Therefore, the patch 
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position (including patch length) is visualized by means of a heat-map, which 
weights the resulting fitness values. The results along the Pareto front are 
summarized in a weighted form, thereby the weighting factor is calculated by 


Wj 
1 . 
& ans) i á ia) i re 
=1..n 


with n as the number of individuals on the Pareto front, Fy max is the maxi- 
mum value observed for fitness 1, and Fy max the maximum value observed 
for fitness 2, respectively. The weighting is normalized and plotted on a uni- 
form mesh, which is only used for visualization purposes, cf. Figure 6-7. 
Thereby, areas with high patch concentration, cf. red regions in Figure 6-7, 
should be preferred over those with a lower patch concentration, cf. blue re- 
gions in Figure 6-7. 


The resulting heat map taken from the application example, cf. Section 7.2, is 
used here to demonstrated the function of the heat map. Therefore, three 
different approaches to create the heat map are compared. The first and most 
evident way, is to visualize only the individuals from the final Pareto front, cf. 
Figure 6-7. With this approach, those areas are visualized, where the patches 
have the largest impact for the reinforcement criteria. 
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Figure 6-7: Heat map with the individuals from the final Pareto front only, where the colour 
represents the patch concentration along the final front (red = high patch concen- 
tration, blue = low patch concentration), while the abscissa and ordinate axes rep- 
resent the initial components x and y coordinates 


The second way is to visualize all individuals incorporated in the optimization 
process, cf. Figure 6-8. This results in a more unclear representation, since a 
significantly larger number of solutions, spread over the whole search space, 
is incorporated in the visualization process. Nevertheless, the advantage of 
this method is the visualization of the areas being searched by the optimiza- 
tion. 


The last visualization method is shown in Figure 6-9, where a user-defined 
number of individuals from the Pareto front is visualized. This approach is 
useful, if the Pareto front can be constrained to a specific range, to be visual- 
ized. With this approach, the visualization results in a clearer picture, since 
only a smaller number of individuals is taken into account. 
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Beside the performance evaluation, the heat map is a useful tool for the de- 
signer to find the areas in which the reinforcement is most effective. While 
the visualization proposed in Figure 6-8 is more useful for the performance 
evaluation, the visualization of a selected number of individuals from a con- 
straint area on the Pareto front is a helpful tool for the designer. 
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Figure 6-8: Heat map with all results obtained during the optimization (colour and axes as for 
Figure 6-7) 
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Figure 6-9: (b) Heat map with a selected number of individuals (colour and axes as for Figure 
6-7); here the 15 stiffest individuals are visualized ((a) red dots in the Pareto front) 
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7 Application examples and numerical 
results 


7.1. Plate structure 


In this section, three plate examples are used to demonstrated the influence 
of the crossover type (plate-1), the influence of different patch numbers 
(plate-2) on the performance of the algorithm, and the behavior of the algo- 
rithm with a large number of generations (plate-3). 


Application example plate-1 


The first application example is used to demonstrate the ability of the method 
to solve the patch optimization problem. Based on equation (9), the optimi- 
zation problem is formulated as 


Minimize F(u) (64) 
= Minimize [global strain energy, patch usage]', 


with the definition of the global strain energy given in equation (58). Further- 
more, the influence of the crossover type is studied with this example. There- 
fore, the setup of the test calculations is shown in Table 7-1. 


The load case is a bending load, cf. Figure 7-1 (b), where the design area is the 
top side of the plate. For the crossover comparison, the basis setup A-Opt 1 
is used as benchmark. To characterize the distribution within the search 
space, all solutions are arranged in groups with patch-length steps of 25 mm, 
from 0 mm to 425 mm maximum patch length for example A. The division 
into the seventeen groups will also be used for A-Opt-2. The distribution of 
the results within the search space is shown in Figure 7-1 (a). The presented 
figure shows that solutions with smaller total patch length are more likely to 
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appear. This is also demonstrated by the percentage share of all A-Opt 1 re- 
sults in Figure 7-1 (d). In the given example, the shortest possible patch length 
is 25 mm, therefore patch usage group 1 remains empty. Furthermore, four 
reference points are given in Figure 7-1 (a), which correspond to the patch 
configurations shown in Figure 7-1 (c). These points shall help to evaluate the 
optimization results. 


Table 7-1: Setup of the evolutionary algorithm used for the test calculations with the plate 


structure 
Name Crossover Muta- Replace- Nich- Number of 
tion ment ing patches 
A-Opt1 Parameterized High Elitist No 2 
binary 
A-Opt2 Multi-point bi- High Elitist No 2 
nary 


Comparing the Pareto fronts, both optimizations converge to a similar front, 
cf. Figure 7-2 (a), which allows the assumption that the proposed approach 
runs stable. The convergence rate, represented by the increase of the domi- 
nated area (Figure 7-2 (b)), is significantly higher for the multi-point binary 
crossover (A-Opt 2) than for the parameterized binary crossover (A-Opt 1). 
This is due to the fact, that the multi-point binary crossover can change a de- 
sign variable on several bits, while other design variables remain unchanged. 
Consequently, the multi-point binary crossover has a more pronounced 
global search than the parametrized binary crossover, which only applies one 
change per variable. Comparing the number of solutions found on the final 
Pareto front, the multi-point binary crossover also has a better performance. 
Both heat maps show the similar result that a reinforcement on the left side 
of the plate will create the best result in terms of the total strain energy, Fig- 
ure 7-2 (e) and (f). This can also be seen by the proposed reference points, cf. 
Figure 7-1 (c), where configuration 4 is the best solution, in terms of strain 
energy minimization. Comparing the distribution of the solutions on the Pa- 
reto front of A-Opt 1 (cf. Figure 7-2 (d)) with the distribution of all solutions 
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of A-Opt 1 (Figure 7-1 (d)), it can be seen that the qualitative distribution is 
very similar. As for A-Opt 1, the solutions along the Pareto front of A-Opt 2 
also show a tendency to smaller total patch lengths. 
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Figure 7-1: Distribution of all calculated individuals in the search space of A-Opt 1 (para- 
metrized binary crossover) in comparison with four reference solutions (a). The blue 
points represent the final Pareto front. Load case used for application example 
plate-1, where A represents the constraint and B the load (b). Patch configurations 
of four reference solutions (c). The red bars illustrate the resulting percentage share 
of all individuals (blue and grey dots in (a)) for each patch-usage group with a step- 
length of 25mm (d). (from [103], modified) 


Summarizing the crossover effects, the multi-point binary crossover offers a 
better performance in terms of convergence rate and number of found solu- 
tions. Since both final Pareto fronts are quite similar, the multi-point binary 
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crossover is used to perform the optimization runs with the curved structure, 
cf. Section 7.2. 
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Figure 7-2: 


Comparison of the results of the two crossover processes A-Opt 1 (parameterized 
binary crossover) and A-Opt 2 (multi-point binary crossover): (a) Comparison of the 
final Pareto fronts, (b) the dominated area for the evaluation of the convergence 
behavior, (c) development of the number of elements on the current Pareto front 
for each generation, (d) percentage share of the solutions along the Pareto front 
(blue dots in (a)) for each length-group with a step-length of 25mm, (e) heat map 
for A-Opt 1, (f) heat map for A-Opt 2. [103] 
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Application example plate-2 


As second application example a plate structure with a different load case is 
used, to demonstrate the influence of the patch number on the final optimi- 
zation result. Furthermore, the influence of the mutation probability is inves- 
tigated. This study was published by Fengler et al. [131]. Based on equa- 
tion (9), the optimization problem is formulated as 


Minimize F(u) (65) 
= Minimize[displacement, patch usage]”. 


With the first fitness value as the displacement, a local criterion is used in- 
stead of a global criterion. The amount of patch used as reinforcement is used 
as second criteria. The load case is presented in Figure 7-3 (b), where the con- 
strained side is represented by A, while the load is applied in form of a line 
load at B. In addition, the line at B is used as evaluation set for the displace- 
ment criterion. 


To demonstrate the influence of the number of patches on the optimization 
result, the patch number is set to Np = [1,2,3], while all other parameters 
remain constant during the optimization run, cf. Table 7-2. The maximum 
number of generations for this example is n = 50. 


Table 7-2: Setup of the evolutionary algorithm used for the patch number influence study 


Name Crossover Muta- Replace- Nich- Number of 
tion ment ing patches 

A-Opt3 Parameterized High Below- No 1 
binary limit 

A-Opt4 Parameterized High Below- No 2 
binary limit 

A-Opt5 Parameterized High Below- No 3 
binary limit 
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The mutation parameter is set to Vm = 0.2 for every calculation, cf. Section 
3.3.2. Based on this and with equation (21) the actual probability for mutation 
can be calculated, depending on the number of patches and thus the number 
of design variables. The calculated actual probabilities for three possible 
events no mutation, mutation of 50 % of the design variables and mutation 
of all design variables are given in Table 7-3. Comparing the probabilities, it 
can be seen, that for the one-patch optimization the majority of individuals 
remains un-mutated, while for the optimization with larger patch numbers, 
the amount of unchanged solutions is considerably lower. A mutation of all 
design variables would comply with a random search, and is therefore unde- 
sirable. Looking at the probability of a mutation in every design variable, it 
can be seen that it is negligible small for the one-patch optimization, and al- 
most zero for A-Opt 4 and A-Opt 5. 


Table 7-3: Calculated actual probability of mutation, based on the number of design variables 
and a probability of mutation Vm = 0.2 


Name No mutation Mutation of 50 % Mutation of all de- 
of the design varia- sign variables 
bles 
A-Opt 3 40.96 % 15.36 % 0.16 % 
A-Opt 4 16.77 % 4.59% ~0% 
A-Opt 5 6.87 % 1.55% ~0% 


The final Pareto sets for the different patch number configurations are shown 
in Figure 7-3 (a). The resulting curves show similarities in fitness ranges that 
could be covered by all patch configurations. These similar results are caused 
by multiple patch solutions, with one long patch and one or two patches with 
the defined minimum patch length. The longest patch is located here in a sim- 
ilar position, as for the calculation with fewer patches. This can also be 
demonstrated with the heat maps, cf. Figure 7-4 (b) to (d), created for the 
section of the Pareto front covered by all patch configurations. 
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Figure 7-3: Load case and evaluation of the example structure plate-2 for different patch num- 


bers: Comparison of the different Pareto fronts (a), load case used for all calculation 
(b), the dominated area for the evaluation of the convergence behavior (c), devel- 
opment of the number of elements on the current Pareto front for each 
generation (d) 


The progress of the optimization run is demonstrated with the dominated 
area, Figure 7-3 (c). Here can be seen, that the dominated area for the solu- 
tions with one patch (A-Opt 3) is, as expected, significantly smaller than for 
the two- and three-patch runs. Furthermore, it can be seen that the differ- 
ence between the two- and three-patch solutions is rather small, which is also 
visualized in Figure 7-3 (a). In addition, A-Opt 3 converges significantly faster 
to a final front, cf. Figure 7-3 (c). This is due to the fact that the number of 
design variables is smaller here. When looking at the number of solutions 
found, the quantity of solutions is the largest for the A-Opt 3. The explanation 
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therefore is, that once the optimization has converged, the focus is on in- 
creasing the number of solutions along the front. 


Comparing the section of patch usage covered by all three optimization runs, 
cf. Figure 7-4 (a), it can be seen that all optimization runs create a Pareto front 
with distinct steps. This is caused by the discrete modeling of the design pa- 
rameter length. When using more than one patch, the gaps between the 
steps are closed, since the total length is a combination of the individual de- 
sign parameter. The values in-between two steps is thereby created by 
patches cut at the border of the plate, which is more likely to appear for larger 
patch numbers. The effect of a significantly larger amount of solutions be- 
comes clear, when comparing the heat maps. Here, A-Opt 3 results in a more 
uniform visualization than A-Opt 4 and A-Opt 5, cf. Figure 7-4. Nevertheless, 
all optimization runs have a similar result. 


Summarizing, the findings from the patch number study, the presented ex- 
ample showed that different patch numbers create similar results, in areas of 
the Pareto front accessible with one, two or three patches, cf. Figure 7-4 (a). 
Furthermore, it could be seen, that the influence of the discrete modeling of 
design parameters is less pronounced for larger patch numbers, since the 
gaps between two steps in the Pareto front are closed by solutions cut at the 
border of the plate, as the Pareto front shows, cf. Figure 7-4 (a). Additionally 
the influence of the number of design parameters on the actual mutation 
probability has been discussed. Thereby has been shown, that for larger patch 
numbers the probability for a mutation of all design parameters is signifi- 
cantly low, and therefore it is ensured that the crossover remains the main 
search operator. 


95 


7. Application examples and numerical results 


8 


s 


F2- Patch usage [mm] 


8 


(a) d a AOR 3 (b) i 

\ —+— A-Opt 4 os 

SET eA Opt 5 os 

160 a 0.7 
te 

4 . os 

: 04 

5 3 o3 

02 

0.4 

o : 7 k 


20 22 24 26 -150 -100 -50 o so 100 150 
Fitness 1 - Displacement [mm] 


1 1 
(C) sao bs (d) bs 
08 os 
or o7 
2 os 06 
05 os 
04 04 
03 o3 
02 loz 
x 01 2 o1 
o o 


-150 -100 50 o so 100 150 -150 -100 50 0 so 100 150 


Figure 7-4: Considered section of the Pareto front (a), and the corresponding heat maps: A- 
Opt 3 (b), A-Opt 4 (c), A-Opt 5 (d) 


Application example plate-3 


For the last plate example, the same load case and the same optimization 
objective as for plate-2 is used, cf. Figure 7-3 (b) and equation (64). Compared 
to plate-2, the maximum number of generations has been raised from Nmax = 
50 to Nmax = 200. However, the optimization converged after nmax = 136 
generations. The further settings of the algorithm are the same as for test 
case A-Opt 5, see Table 7-4. 


Table 7-4: Setup of the evolutionary algorithm used for the study of the influence of the max- 
imum number of generations 


Name Crossover Muta- Replace- Nich- Number of 
tion ment ing patches 
A-Opt6 Parameterized High Below- No 3 
binary limit 
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The final Pareto front found during the optimization is compared to A-Opt 5 
to evaluate its quality, Figure 8-4 (a). It can be seen, that the optimization with 
a higher number of generations can increase the performance of the Pareto 
front, especially in the area of higher patch usage. This is also visualized by a 
significantly larger dominated area, cf. Figure 8-4 (b). Since the largest part of 
the performance increase is within an early generation, it is obviously not 
caused by the increase of the generation number. The explanation for the 
increase is the mutation, which is a random based procedure and can there- 
fore randomly vary between a local search, as for A-Opt 5, and a global 
search, as for A-Opt 6, where the global search is characterized by larger im- 
provements within the dominated area between two successive generations. 


As expected, the number of found solutions for the optimization with the 
higher number of generations, A-Opt 6, is significantly higher than for A- 
Opt 5. However, the number of solutions found up to the 50" generation is 
approximately similar. Comparing the heat maps, cf. Figure 7-4 (d) and Figure 
7-5 (d), can be seen that the difference is negligible. 


Concluding the study on the influence of the maximum number of genera- 
tions, the highest impact of a higher number of generations is on the total 
quantity of solutions found on the final Pareto front. Therefore, the domi- 
nated area could be a useful criteria for the convergence of a solution, since 
for the event of convergence of the optimization, the dominated area approx- 
imates a horizontal line. Comparing the qualitative results, in terms of the 
heat map, the influence is negligible, therefore it should be sufficient to con- 
tinue with a smaller number of maximum generations. 
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Figure 7-5: Comparison of the results for the variation of maximum number of generations (A- 
Opt 5: n=AAA, A-Opt 6: n=200), (a) Comparison of the final Pareto fronts, (b) the 
dominated area for the evaluation of the convergence behavior, (c) development of 
the elements on the current Pareto front for each generation, (d) heat map for the 
final generation of A-Opt 6. 
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7.2. Slightly curved structure 


The second example is a slightly curved structure, in which the kinematic 
draping simulation has a larger effect than in the previous example, where 
the kinematic algorithm only checked a possible overlapping of patches be- 
yond the plate boundaries. Likewise as for the first example, the optimization 
problem can be formulated based on equation (9) as 


Minimize F(u) = (66) 
Minimize [global strain energy, patch usage]. 


Based on the findings from the plate structure, the crossover type is kept con- 
stant for this example (multi-point binary crossover), while the optimization 
parameters for mutation, replacement and niching are systematically varied 
to show their influence on the performance of the optimization. The settings 
for the calculation are summarized in Table 7-5. The following results have 
been published by Fengler et al. [103]. 


Table 7-5: Setup of the evolutionary algorithm for the second application example 


Name Crossover Muta- Replace- Nich- Number of 
tion ment ing patches 

B-Opt 1 Multi-point bi- Low Elitist No 2 
nary 

B-Opt 2 Multi-point bi- High Elitist No 2 
nary 

B-Opt 3 Multi-point bi- Low Below- No 2 
nary limit 

B-Opt 4 Multi-point bi- Low Elitist Yes 2 
nary 


The load case is a bending load with two constrained areas, cf. Figure 7-6 (c). 
As for the plate structure, the Pareto fronts obtained for the different setups 
(cf. Table 7-5) are compared with predefined reference points to evaluate the 
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results, cf. Figure 7-6 (a). Thereby, the patch reinforcement of reference point 
4 ends at the constraint area, while the reinforcement from point 3 is limited 
by the part geometry, Figure 7-6 (b). 


Comparing the different mutation settings, the higher mutation (B-Opt 2) cre- 
ates a slightly better result regarding stiffness maximization for large patch 
lengths, cf. Figure 7-7 (a). This is due to the fact that the changes applied by 
the more flexible mutation operator work similar to the more flexible crosso- 
ver process, and therefore increase the search power of the optimization. The 
dominated area does not converge faster compared to the low mutation (Fig- 
ure 7-7 (b)), but starts at a higher value. When looking at the number of found 
solutions, as a second quality criterion, the number of found solutions at the 
final front is higher for the higher mutation rate, cf. Figure 7-7 (c). Comparing 
the heat maps, it can be seen, that both optimization runs end up in similar 
results, cf. Figure 7-10 (a) and (b). The better performance in terms of finding 
the extreme point regarding stiffness maximization is not visible in the heat 
map (Figure 7-10 (b)), since only a small proportion of solutions is in the ex- 
treme area, and is therefore underrepresented within the heat map. This 
shows also, that the heat map is an additional tool, but should not be consid- 
ered as the only basis of decision, but always in combination with the Pareto 
front. 


Comparing the results obtained by calculation B-Opt 1, B-Opt 3 (replacement 
type below-limit, where also dominated individuals are passed to the next 
generation), and B-Opt 4 (with niching), no significant change regarding the 
Pareto front and the convergence rate can be seen, cf. Figure 7-8 and Figure 
7-9 each (a) and (b). The difference in the results of the different setups be- 
comes more obvious comparing the heat maps, cf. Figure 7-10. Here it can be 
seen, that the distribution of the resulting patch positions is more balanced 
(reflected in a more pronounced symmetrical picture in the heat map), when 
more fronts are taken into account (B-Opt 3) or a niching pressure is applied 
(B-Opt 4). This effect results from a higher diversity of the found solutions 
along the Pareto front. 
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Figure 7-6: (a) Comparison of the Pareto fronts obtained for the curved structure, (b) prede- 
fined reference points to classify the results, and (c) load conditions. (from [103], 
modified) 


For the given example, a reinforcement covering the position of the structural 
kink has a large effect on the total component's stiffness. This becomes evi- 
dent when comparing reference points 2 and 3, cf. Figure 7-6 (b), as well as in 
the heat maps, represented by the red areas, cf. Figure 7-10 (c) and (d). The 
red areas in the heat maps for B-Opt 1 and B-Opt 2 are smaller compared to 
those of B-Opt 3 and B-Opt 4, resulting from a lower diversity, and therefore 
an allocation of solutions with smaller patch length. 
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Summarizing the results obtained with the curved structure, all optimizations 
converge to a similar final Pareto front. This shows that the proposed ap- 
proach is capable to solve the patch optimization problem. Furthermore, it 
could be seen that a higher mutation rate is likely to improve the ability of 
finding extreme points regarding stiffness maximization. However, the diver- 
sity of the solutions is still low, which results in a more asymmetrical heat 
map. Utilizing the heat maps, the best results could be achieved by incorpo- 
rating a niching pressure and using the below-limit replacement method. This 
was demonstrated by a more symmetrical picture in the heat map. The effect 
on the convergence rate was found to be small for the proposed configura- 
tions, except for the higher mutation rate, which creates a slightly higher con- 
vergence rate, represented by the dominated area. Based on this, a higher 
mutation rate is suggested, if the finding of extreme solutions (maximum stiff- 
ness) has a higher preference. Additionally, a niching pressure and a below- 
limit replacement are helpful to achieve a more balanced Pareto front and 
therefore more suitable heat maps. 
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Figure 7-7: Comparison of the results for the variation of the mutation (B-Opt 1: low mutation 
rate, B-Opt 2: high mutation rate), (a) Comparison of the final Pareto fronts, (b) the 
dominated area for the evaluation of the convergence behavior, (c) development of 
the elements on the current Pareto front for each generation, (d) percentage share 
of the solutions along the Pareto front for each length-group [103, 103] 
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Figure 7-8: Comparison of the results for the variation of the replacement type (B-Opt 1: elitist 
replacement, B-Opt 3: below-limit replacement) , (a) to (d) as inFigure 7-7. [103] 
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Figure 7-9: Comparison of the results for no-niching and niching (B-Opt 1: no-niching, B-Opt 4: 
niching), (a) to (d) as in Figure 7-7. [103] 
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Figure 7-10: Top view of the generated heat-maps for the results obtained along the Pareto front 


4 


(red: high patch concentration, blue: low patch concentration): (a) basic setup (B- 
Opt 1), (b) higher mutation rate (B-Opt 2), (c) change of replacement type (B-Opt 
3), (d) niching pressure (B-Opt 4). [103] 


106 


7.3. IRTG reference structure 


7.3. IRTG reference structure 


In this section the reference structure developed within the DFG-GRK2078 
will be used to demonstrate the performance of the algorithm. First a stiffness 
optimization is performed with both crossover operators, to demonstrate the 
influence of the crossover for a rather complex structure (example 1). The 
second optimization run consists of a combined warpage and stiffness opti- 
mization along the proposed optimization chain (example 2), cf. Section 6.2, 
for which only one crossover operator is considered. In the presented exam- 
ples, the fiber orientation from a mold filling simulation is used. 


Application example GRK2078-reference structure 1 


Based on equation (9) the optimization problem is formulated as 


Minimize F(u) = (67) 
Minimize[strain energy, patch usage]. 


The parameters used for the evolutionary algorithm are given in Table 7-6. 
Here the crossover operator is varied to validate the findings from the plate 
structure, cf. Section 7.1, with a 3D structure. Furthermore, a niching pressure 
is applied here, since this is excepted to result in a more uniform distribution. 
In both cases the optimization is run for n = 30 generations. 


Table 7-6: Setup of the evolutionary algorithm for the study with the 3D structure 


Name Crossover Muta- Replace- Nich- Number of 
tion ment ing patches 
C-Opt1 Parameterized Low Elitist Yes 2 
binary 
C-Opt2  Multi-point bi- Low Elitist Yes 2 
nary 
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The load case used for the structural simulation is a four-point-bending load, 
cf. Figure 7-11 (a), where the structure is loaded with two line loads at B and 
constrained at four support points Al to A4. Additionally, it should be men- 
tioned that the proposed structure is not symmetric. The slope of the two 
sides C and D, cf. Figure 7-11 (b), is 140° and 120°. Whereas the two sides E 
have an equal slope of 135°. 


\ —_ 


Figure 7-11: Reference structure used as a 3D validation structure: (a) load case with the four 
support points A1 to A4 and the line loads B, (b) top view of the structure with the 
sides with different slope C and D, and the sides with an equal slope E 


Comparing the Pareto fronts, the multi-point crossover creates a slightly bet- 
ter result than the parametrized binary crossover. This has also been ob- 
served in the first application example, cf. Section 7.1. Furthermore, the un- 
symmetrical behavior of the example structure can be seen by looking at the 
reference points. When comparing reference point 3 and reference point 4, 
cf. Figure 7-12, the difference in the performance of mirrored solutions can 
be seen. This effect is enlarged by the process-induced fiber orientations, 
which are due to the asymmetric geometry also asymmetric. 
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Figure 7-12: Comparison of the Pareto fronts achieved with the two crossover types (C-Opt 1: 
parameterized binary, B-Opt 2: multi-point binary) with four reference points to 
classify the results. In reference point 3 and 4 the smaller patch is equal to the min- 
imum patch length of 25 mm. 


The performance of the two crossover types is compared with the develop- 
ment of the dominated area, cf. Figure 7-13 (a). It is noticeable, that unlike in 
the previous examples, a large increase of the dominated area between gen- 
eration 1 and 2 occurs. Additionally, for both crossover types, there are 
phases in which non or almost non progress is achieved, like between gener- 
ation 2 to 10 for C-Opt 1 or generation 6 to 14 for C-Opt 2. This can be ex- 
plained by the probability based characteristic of the evolutionary algorithm, 
which can results in phases with large- or non-progress, as in the presented 
example. Comparing the number of individuals on the front per generation, 
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cf. Figure 7-14 (b), C-Opt 2 finds also a larger number of solutions. However, 
C-Opt 1 creates a significantly larger number of solutions in the group of 


larger patch lengths. Compared to the previous examples, cf. Section 7.2, the 
distribution within the search space is more balanced here. The influence of 
the fiber orientation and of the asymmetrical geometry is visible in the heat 


maps, cf. Figure 8-4. In both cases, the area with a more pronounced rein- 


forcement is on the side with the smaller angle. 
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binary, B-Opt 2: multi-point binary), (a) Comparison of the dominated area for the 
evaluation of the convergence behavior, (b) development of the elements on the 
current Pareto front for each generation, (c) percentage share of the solutions along 
the Pareto front for each length-group 


7.3. IRTG reference structure 


Summarizing the results for the optimization of the GRK2078 reference struc- 
ture, both crossover types achieve a similar result in terms of the final Pareto 
front, which shows that the proposed algorithm is able to solve the patch op- 
timization problem for a complex 3D structure independent of the evolution- 
ary configuration. Again, the performance of the multi-point crossover is bet- 
ter than the one of the parametrized binary crossover, which was also 
observed in the first example. 


(a) (b) ” 


(c) (d) ” 


Figure 7-14: Heat maps of the two optimization runs, each with perspective and top view: (a) 
and (b) with parameterized binary crossover (C-Opt 1); (c) and (d) with multi-point 
binary crossover (C-Opt 2); 
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Application example GRK2078-reference structure 2 


The last application example is used to demonstrate the performance of the 
proposed algorithm, by incorporating the warpage simulation in the optimi- 
zation loop. For this purpose, the same geometry as for the previous example 
is used. Utilizing equation (9), the optimization problem can be formulated as 


Minimize F(u) (68) 
_ Min [strain energy (structural load case), 
~ strain energy (curing), distortion]" 


Hereafter, the identifiers F}, F2, and F} are used for strain energy resulting 
from the structural load case, strain energy resulting from curing, and the dis- 
tortion objective. For the structural load case, the same conditions are ap- 
plied as for reference structure 1, cf. Figure 7-11 (a). The optimization of the 
occurring warpage is performed with two objectives, minimization of process- 
induced displacements and strain energy. The reason for considering both is 
that a component can be distortion-free, with the drawback of a significant 
amount of residual stresses. For the distortion minimization objective, the 
mean value from the four evaluation points C1 to C4 is used, Figure 7-15 (a), 
while the strain energy is taken from the whole model. The temperature pro- 
file used for the curing and warpage simulation is presented in Figure 7-15 (b). 


The setup of the evolutionary algorithm is the same as for C-Opt 2, and is 
shown in Table 7-7. 


Table 7-7: Setup of the evolutionary algorithm used for the warpage optimization 


Name Crossover Muta- Replace- Nich- Number of 
tion ment ing patches 
C-Opt 3 Multi-point bi- Low Elitist Yes 2 
nary 
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Figure 7-15: Load case for the warpage simulation: (a) model with the constraint points (A1: 
x and y direction, A2: z direction) and the evaluation points for the distortion crite- 
ria (C1 to C4), (b) temperature profile used for the curing and warpage simulation 


The three-dimensional Pareto front, as a result of the combined structural 
and warpage optimization, is shown in Figure 7-16 (a). In addition to the Pa- 
reto optimal results, an approximation of the Pareto front is shown. Besides, 
the Pareto front is analyzed in each search space individually, cf. Figure 7-16 
(b) to (d). Therefore, all results from the 3D Pareto front are shown in the 
individual search spaces, and the subset which would correspond to the 
search space Pareto front is highlighted. 


The most distinct creation of a front is shown for the search space covering 
the two strain energy objective functions, cf. Figure 7-16 (b). This is due to the 
fact that the two objective functions are most contradictory, since longer 
patches are required for higher stiffness, which also leads to higher curing 
stresses. Looking at the distribution in the F2/F3 search space, it can be seen 
that the majority of solution is allocated in one area the search space. This is 
caused by the close link between the two objective functions. Nevertheless, 
this is not the fact for all solutions, since there are solutions for which both 
objective functions differ considerably. These results correspond with solu- 
tions with ether high distortion and therefore low residual stresses or vice 
versa, therefore it is reasonable to consider both objective functions during 
the optimization, to ensure a solution with a minimum of distortion and re- 
sidual stress. 
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Furthermore, the result of the optimization is visualized with a heat map, cf. 
Figure 7-17. Comparing the heat map with those from the first application 
example, cf. Figure 7-14, can be seen that the highest patch concentration is 
on the same side, though the area with a high patch concentration is smaller 
when considering the warpage objective during the optimization. Besides 
this, the patch length is in general small than for the previous example. This 
is caused by the two manufacturing process based objectives, F2 and F3, since 
a smaller patch length results in lower distortions and residual stresses. In 
addition, the patch concentration on the side with the higher geometrical 
stiffness, cf. Figure 7-17 (b) lower part of the component, is higher. This leads 
to the deduction that the ratio of stiffness improvement to resulting warpage 
inthis area of the component is better than the ratio of stiffness improvement 
to patch usage, and thus the use of patch at this point is more advantageous, 
with regard to the optimization objectives for the combined warpage and 
stiffness optimization. 


Summarizing the findings of the warpage optimization, a three-dimensional 
Pareto front could be determined, considering strain energy from a structural 
load case and, distortion and strain energy from the curing process as objec- 
tive functions. Based on the given example, the combination of the three ob- 
jectives is a useful approach for the optimization of local reinforcements, 
since manufacturing influences can be considered during the optimization 
process. In addition, the heat maps are again a very good way to visualize the 
optimization results. 


114 


7.3. IRTG reference structure 


(a) 


25 


3 
8 
4 
8 
g 
& 
5 
5 
5 


4 
F2 - Strain energy (curing) {J} 6 7 
3 4 z 
0 ‘ 2 
F3 - Distortion [mm] 
x 8 
(b) <] Ham (0 o] 
4 |e Pareto front F1/F2 1 = 
= 504 | 5 
Sul e+ 
ə = 
* A 
5 & | . 
3 A . . 
5 54 
3 C s ‘7 
$ rd 
£ t 4 6 | 
j t . e | 
@ 25 . s u, 
«4 4 Ki, Aa 7 
£204 | Aa. | 
154 | 
| $ i H I H 5 o4 r - 1 - + ; 
15 20 25 30 15 20 25 30 
F1 - Strain energy (structural load case) [J] F1 - Strain energy (structural load case) [J] 
(d) » (e) 
. © Pareto Font (al fines values) 042 set 
e- Paret font F21F2 RER Per ezi 
_s “| | | | | et” 
B » Fow / 
3 : . 8 r 
< . 
54 3 s 
2 d| o gos f 
f ATTS 5 / 
z ~? P 
2 f 
2 I 0.36 } 
P: 
« 
o + om 
15 20 25 30 35 40 45 50 55 o 5 0 5 2 5 w 
F2- Strain energy (curing) [J] Generation [-] 
100 et 
Met 
.. 
z% y 
E 
È a Z 
GJ IE 
5 “ k 
g pae. 
È 0 
w « 
. 
2 
o t t t t 
o 5 10 5 a 2 30 
Generation [-] 


Figure 7-16: 


Results for the combined structural and warpage optimization: (a) 3D representa- 
tion of the final Pareto front (results and interpolation); Analysis of the Pareto front 
in the three solution spaces (b) strain energy from curing and load case, (c) strain 
energy from load case and distortion, (d) distortion and strain energy from curing, 
(e) dominated area for the evaluation of the convergence behavior, (f) development 
of the elements on the current Pareto front for each generation 
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(a) (b) "| 


Figure 7-17: Heat maps of the optimization run with warpage consideration, (a) perspective 
and (b) top view 


7.4. Conclusion on application results 


Three application examples have been studied to show the ability of the pro- 
posed algorithm to solve the patch optimization problem. Starting from a ra- 
ther simple plate structure, the geometrical complexity was increase to a 
slightly curved and finally to a 3D reference structure. 


The plate structure was used to perform basic feasibility studies and to com- 
pare the performance of the different crossover operators. Thereby, the 
multi-point binary crossover offers a better performance in terms of conver- 
gence rate and number of found solutions, while the found Pareto fronts are 
similar with both crossover types. Furthermore, the influence of the patch 
number on the performance and quality of the found Pareto fronts has been 
analyzed. Thereby, it could be seen that gaps in the Pareto front, caused by 
the discrete modeling of the design parameters, are closed when using a 
larger n umber of patches. The reason therefore are patches cut atthe border 
of the plate, which are more likely to appear for larger patch numbers. In ad- 
dition, the different patch numbers create similar results, in areas accessible 
with one, two, and three patches. The last study with the plate structure an- 
alyzed the influence of the maximum number of generations on the optimi- 
zations performance. Thereby, could be demonstrated that once the Pareto 
front has converged to a final Pareto front, the highest impact of a higher 
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number of generations is on the total quantity of solutions found. Therefore 
it should be sufficient to continue with a smaller number of maximum gener- 
ations, though it has to be ensured that the optimization converges. 


The curved structure was used to analyse the influence of different settings 
of the evolutionary algorithm on the optimizations performance. Thereby, all 
optimizations converge to a similar final Pareto front, which shows that the 
proposed approach is capable to solve the patch optimization problem. Addi- 
tionally, it could be demonstrated that a higher mutation rate is likely to im- 
prove the algorithm’s ability to find solutions in extreme areas, like high stiff- 
ness or low patch usage. A more balanced distribution of the solutions along 
the final Pareto front can be ensured when a niching pressure is applied. Fur- 
thermore, the influence of the two different selection types, elitist and below 
limit, on the optimization result has been studied. Thereby, the below limit 
method results in a more balanced Pareto front than the elitist method. The 
effect on the convergence rate was found to be small for the proposed con- 
figurations, except for the higher mutation rate, which creates a slightly 
higher convergence rate. 


The last application example studied here, was the reference structure devel- 
oped within the DFG-GRK2078. Again, the performance of the multi-point bi- 
nary crossover was better than those of the parametrized binary crossover. 
After the structural optimization was demonstrated successfully with the al- 
gorithm, a combined warpage and stiffness optimization was performed. 
Contrary to the first example, the three objective functions, strain energy 
from a structural load case, distortion, and strain energy from curing, have 
been used here. Thereby could be demonstrated that the concurrent use of 
distortion and strain energy from curing is a good approach to consider ef- 
fects resulting from warpage in the optimization. 
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8.1. Robustness measurements 


According to Jin and Branke [94], a robust design should still work satisfactory 
when inaccuracies, e.g. from the manufacturing process, influence the design 
parameters. Therefore, perturbations and changes are applied on each design 
variable of the optimal solutions, once the optimization procedure is done. 
This definition of a robust solution is also used in the context of this thesis. 
Therefore, the robustness evaluation is used to rate the solutions on the Pa- 
reto front in terms of sensibility against inaccuracies of the design parameter, 
while the determination of the robust front is not focused here. Methods to 
determine the final robust front have been presented by Deb and Gupta 
[132], and Gaspar-Cunha et al. [93]. 


The variation of the design parameter is applied in the design space, while the 
evaluation of the robustness is done in the solution space, cf. Figure 3-2. 
Thereby, the same variation of the design parameter can result in different 
variations in the solution space, cf. Figure 8-1. Furthermore, small changes in 
the design space can result in larger changes in the solution space. The inac- 
curacies of the design parameters are taken from the manufacturing process, 
cf. Section 6.1, while the allowed deviation in the solution space comes from 
tolerance management and is a user defined input parameter. Based on this 
input, two different ratios are applied to rate the robustness of a solution, the 
degree of robustness and the robustness index. 
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Design space Solution space 


Figure 8-1: Effect of the variation of the design parameters on the solutions for two exemplary 
designs, where equal variation of the design parameters may result in different var- 
iations in the solution space 


8.1.1. Degree of robustness 


The concept of a degree of robustness for multi-objective optimization was 
first introduced by Barrico and Antunes [133]. They defined the degree of ro- 
bustness of a solution u as a value k, for which 


1. the percentage of solutions in the neighborhood kô is greater or 
equal than a predefined value p, and 

2. the percentage of solutions in the neighborhood (k + 1)Ö is lower 
than p. 


Thereby, ô represents the initial test neighborhood size, and kô the size of 
the kt” neighborhood, for which a predefined percentage of solutions Pk, 


ne Ak robust 


Pk , (69) 


Ax sample 


has to be within the neighborhood. Here, Ak, sample is a randomly created set 
of h solutions u; in the kô neighborhood, and Ax robust the set of robust so- 
lutions u, within the same kô neighborhood defined by 


Ax robust (70) 
= {uj € Ax sampie|Fi(uj) < Fio(Auo), Vi} with j 
Dh, 
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with the number of fitness functions i, the test sample points uj, and the al- 
lowed deviation of the fitness function F,,(Au,). In addition, the size of the 
test sample is calculated by 


h=h+(k—-1)*d*ho, (71) 


where hy is the user-defined initial sample size, and d the grow rate for the 
sample size for increasing degree of robustness k. The definition of the de- 
gree of robustness in the design and solution space, cf. Section 3.2, is visual- 
ized in Figure 8-2. 


Figure 8-2: Definition of the degree of robustness k in the design (left) and solution space 
(right) where AF9, and AF 92 represents the maximum allowed deviation in the so- 
lution space ([133] modified) 


8.1.2. Robustness index 


The second robustness measure used here, is the robustness index, and was 
introduced by Gunawan and Azarm [134] for single-objective optimization. 
The method was later extended for multi-objective optimization by Li et al. 
[135], and Gunawan and Azarm [136]. The method applied in this work is a 
variation, proposed by Li et al. [137], where they use the variation of the de- 
sign parameter for the calculation of a sensitivity region. 
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The basis for the robustness index is the calculation of the sensitivity region 
of a test solution ug. The sensitivity set Q is given by 


A(uy + Au) = (72) 
{Au E R|[F (u + Au) — F (u0)]? < [AFy(Auy)]*}, 


where F (uo) is the fitness of the initial point, F (uo + Au) is the current test 
point’s fitness, Au the deviation of the design variables, and AFọ (Auo) the 
allowed deviation of the fitness function, which is called acceptable objective 
variation range (AOVR). 


A consideration of a worst case scenario regarding the sensitivity region is 
assumed, to calculate the worst-case-objective-sensitivity-region (WCOSR), 
cf. Figure 8-3. Here, a worst case scenario means, that always the most sensi- 
tive direction is considered for the calculation of the robustness index. The 
WCOSR method will thereby always result in an underestimation of the actual 
sensitivity region. 


The tolerance region is a hyper-rectangle in the design space, formed by all 
acceptable values of Au around the current point u,, which is defined by 
Aulower < Au < Au“? In Figure 8-3 the tolerance region is visualized ex- 
emplarily for a design space with two design parameters. 


F, AOVR AFoo 


Aus 


Tolerance Region 
upper 
Aus 


Ayre? 
Lẹ Au: 


Figure 8-3: Definition of a tolerance region in the design space (left) and the corresponding ob- 
jective sensitivity region (OSR) with the resulting worst-case-objective-sensitivity- 
region (WCOSR) and the acceptable objective variation range (AOVR) (right), ([137] 
modified) 
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The most sensitive direction, used for the WCOSR, can be found by solving 
the following optimization problem 


1/2 


max Rı(Au) = arcane (73) 
i=1 


with the number of fitness functions n, the radius of the sensitivity region R, 


and 
Fio + Au) — F;(u 
AF;(Au) = Fi(up + Au) — F; (uo) (74) 
AF; o 
Aylower < Au < Awuprer, (75) 


Based on the radius if the OSR, the robustness index can be defined as [137] 


Rr 


n 


where R; is the external radius of the required size of the AOVR. A solution is 
thereby defined as robust, by the mean of the robustness index, if its OSR is 
completely enclosed by the AVOR, cf. Figure 8-3. 


8.2. Kriging meta-model 


The basic definition of ameta-model is given in Section 2.8. Based on this, the 
Kriging meta-model approach consists of a known set of approximation func- 
tions g and a realization of a stochastic process Z. With this, the unknown 
function of interest can be defined by 


fu) = g (wp +Z(u), (77) 
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where beta represents the vector of the estimated trend function coeffi- 
cients. Typical approximation functions, used for the Kriging method, are lin- 
ear, polynomial or constant. A universal Kriging approach can consist of a set 
of approximation functions. Within this thesis, linear, quadratic, and reduced 
quadratic approximation functions are used. [138] 


The covariance matrix of Z(u) is defined by 
Cov[Z(u'), Z(w)] = o?R([R(u‘, w)]) (78) 


with the variance o”, the correlation function Rfor any of the n, sample 
points u! and u), and the symmetric correlation matrix R. Thereby, R is a 
(n, X n,) symmetric matrix with ones on the diagonal. The correlation func- 
tion R is a user defined function. In the presented approach, a Gaussian cor- 
relation function of the form 


Nav 


iyi) — ha i_ iy? 
R(u‘,u/) = exp (ui — ul) * Oy (79) 
k=1 


is used. Here, 9, are the unknown correlation parameters, used to fit the 
model, na, is the number of design variables, while ul and u) are the kth 
components of the sample points u! and u. 


The predicted estimates F for a point î, other than the sample points u, is 
given by 


F(@) = g™(w)B + r!(MR'(f- GB) (80) 


where f is a vector of the size n,, consisting of the sample points, and G is the 
matrix of regression functions. Furthermore, r is the correlation vector be- 
tween the untried point î and the sample points u! , and is given by 


r'(@) = [RG, ut), R(@, u?), ..., RG, w"s)]. (81) 
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In addition, B can be estimated by 


B = (G™R-'G)-1G™R-'f. (82) 


The variance of an untried point 6?(@) of the underlying model ß and f can 
be estimated by 


62(@) = o? |1 —-[g™(@ r™(@)] es a leak (83) 


where 0? represents the maximum variance, and is defined by 
aAT__ A 
a? = (t/n )(f - GB) R- (F - GB) i 


The Kriging model used here is implemented in the Surfpack software library 
[139], which is part of the Dakota software package [130]. 


8.3. Implemented workflow 


The method presented here, was developed within the master thesis of Bas- 
tian Schäfer [140]. 


The general workflow of the developed robustness evaluation method is 
shown in Figure 8-4. The results of the optimization process, cf. Chapter 1, 
can be used as input for a subsequent robustness study. Thereby, not only 
the final Pareto front is used, but further solutions created by the evolution- 
ary algorithm should be part of the training data set. Hence, the first step is a 
data treatment to select suitable data points. Thereby the results are sorted 
according to their domination value, cf. Section 3.4. Furthermore, duplicates 
are removed from the data set. Duplicates can arise, since the patch is cut at 
the border of the initial component. Therefore different values for the design 
parameter length l can result in identical fitness values. It should be defined 
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here, that for all individuals with the same fitness value for patch usage, only 
the individual with the shortest initial length is kept. 


The next step is the selection of the training data, used to create the meta- 
model, cf. Section 8.2. Here, only data points close to the Pareto front are 
taken into account, since the robustness is only evaluated for the points on 
the final Pareto front, and therefore the meta-model creation focuses on this 
part of the design space. For this purpose, the acceptable objective variation 
range (AOVR), which corresponds to the range in which the solutions are con- 
sidered to be robust, is multiplied by a scaling factor, to define the range 
which is used as training set. Based on the training data, the meta-model is 
created. Therefore, the data is divided into ten equal sized sets, where nine 
sets are used to train the model, and one set is used for verification, to predict 
the resulting model error. Here, every combination of the ten data sets is used 
to train a model, and the mean error value of all ten models is used as quality 
criteria for the model. This procedure is repeated for different trend functions 
of the Kriging model, cf. Section 8.2. Finally, a meta-model for each objective 
function is selected, and trained with all training points within the defined 
training data set. 


The initial meta-model, created before the training data set was increased, is 
used at first to calculate the robustness index, cf. Section 8.1.2. Afterwards, 
the model construction and calculation of the robustness index is repeated 
for the meta-model of the increased data set. This step is necessary to check, 
if the number of training points for the meta-model was sufficient enough. 
Hence, the quality of the model is rated by comparing the calculated robust- 
ness index changes, resulting from the increased training data set, with a pre- 
defined threshold. Since the calculation of the robustness index is much faster 
than for the degree of robustness, the robustness index is used here. After 
the meta-model has been setup successfully, the degree of robustness can be 
calculated as well, cf. Section 8.1.1. 


After the calculation of degree of robustness and robustness index, the results 
are visualized. Besides a diagram for each robustness measurement, a heat 
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map is created. The heat map is created similar to those presented in Sec- 
tion 6.4. Here, the degree of robustness is mapped on the structure, thereby, 
only one solution is considered at a time. Besides the original solution, the 
areas for different values of k are mapped on the initial component. 


The presented workflow of the robustness evaluation is implemented in 
MatLab. For the meta-model, an implementation in Dakota [130] is used. 
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Optimization results 


MatLab Dakota Tool Box 


Data Treatment 
Training data selection 


Meta-model training 


Cross-Validation 
Cross-Validation error 


Meta-model selection ET Meta-model training 
Meta-model | LO 
Robustness index 
calculation 


Evaluate iterative 
robustness change 
Change < Threshold 
or max. iteration 
Degree of robustness 
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Visualization 


Robustness index 


== | Ei m 


Figure 8-4: Workflow of the robustness evaluation method, based on the results created by 
the optimization [140] 
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8.4. Application example 


As application example for the robustness evaluation workflow, the results 
from the optimization run B-Opt 4, presented in Section 7.2, are taken. 


The parameter used for the calculation of robustness index and degree of ro- 
bustness are given in Table 8-1 . The allowed deviations of the fitness values 
AF; o and AF; are given as absolute values, instead of a percentage change, 
since this definition is more applicable here. For example for the patch-usage 
fitness criterion, a definition in percent would result in larger acceptable 
length changes for longer patches, which is not the case with an absolute de- 
viation definition. 


Table 8-1: Parameter setup for the robustness evaluation of the application example 


Degree of robustness 


Parameter Symbol Value Unit 
Neighborhood size in x Ox +1 mm 
Neighborhood size in y by +1 mm 
Neighborhood size in So +1 ? 
Neighborhood size in l ô +1 mm 
Initial sample size ho 100 - 
Grow rate d 0.5 - 
Threshold p 95 % 
Robustness index 

Parameter Symbol Value Unit 
Deviation of design variable x Au, +3 mm 
Deviation of design variable y Au, +3 mm 
Deviation of design variable p Aug +3 ? 
Deviation of design variable l Au, +3 mm 
AOVR fitness 1 AF, 9 + 0.004 J 
AOVR fitness 2 AF, o +10 mm 
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The robustness index 7, calculated for each individual on the final front is 
shown in Figure 8-5 (a). In general, solutions with 7 < 1 are considered as 
robust solutions, since here the maximum deviation of the design parame- 
ters, defined in Table 8-1, results in a lower deviation of the fitness values 
than the applied constraint. The boundary between the robust and non-ro- 
bust region is furthermore marked with the black solid line (7 = 1). In addi- 
tion, the degree of robustness was calculated for each individual on the Pa- 
reto front, cf. Figure 8-5 (b). Here a higher k value corresponds to a higher 
robustness of a solution. 


A degree of robustness of k = 3 is comparable with a robustness index of 
7 = 1, considering the presented example parameter set. The slight differ- 
ence is caused by the different approaches to calculate the robustness meas- 
urements. When calculating the degree of robustness, only a defined thresh- 
old p (in this example 95%) has to be within the allowed deviation, while the 
robustness index calculates the worst case with a given deviation, which 
would be equal to a threshold of p = 100%. 


Comparing both robustness measurements, it can be seen that the predicted 
robust regions have a similar trend. In both cases the solutions at the extreme 
ends of the Pareto front, high stiffness or low patch usage, have a lower max- 
imum robustness (green dots in Figure 8-5). This effect could be caused by 
the meta-model, because the quality of the meta-model, used to calculate 
the robustness measurements, is lower in the extreme areas, since a lower 
number of training points in present here. In addition to the robustness of 
each solution, the Pareto front is plotted with black dots, which helps to clas- 
sify the results. Furthermore, the scattering of the degree of robustness, in 
the middle region of the Pareto front, results from solutions for which only 
small deviations in the design parameters cause large deviations of the fitness 
values. 
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(a ) Analysis of the robustness index for the first Pareto front 


patchlength [ - ] 0 0 displacement [- ] 


* projection of the first Pareto front onto the response-plane 


Analysis of the degree of robustness for the first Pareto front 


k[-] 


300° 1 


patchlength [ - ] 0 0 displacement [ - ] 


* projection of the first Pareto front onto the response-plane 


Figure 8-5: Results of the calculation of robustness index (a) and degree of robustness (b), 
both in comparison with the final Pareto front (black dots) 


Based on the degree of robustness, a heat map is created. The visualization 
via a heat map is only meaningful for the degree of robustness, since here a 
distinction between the different robustness neighborhoods k can be made, 
while for the robustness index, only the area which was a-priori defined to be 
robust would be visualized. 
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In Figure 8-6, the point with the best performance regarding stiffness, i.e. fit- 
ness 1 (b) and a point with a high robustness (c) are used for demonstration 
purposes. The area marked in the heat maps has a similar size for both cases. 
However, the number of neighbourhood enlargements, and thus the degree 
of robustness, is significantly higher for point 2 than for point 1, cf. Figure 8-6 
(b) and (c). By comparing the robust areas it can be seen, that changes of the 
patch position affect the fitness values more than changes in the orientation 
of the patch. This is demonstrated by a larger robust area, caused by orienta- 
tion changes. Furthermore, the influence of the chosen draping method can 
be seen in the heat maps, since all changes of the orientation are performed 
with the center of the patch as center of rotation. 
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Figure 8-6: (a) Pareto front with the selected points and the created robustness heat maps, (b) 
for the stiffest point (point 1, k = 4), and (c) for a point with a high robustness 
(point 2, k = 10) 


132 


8.5. Conclusion on robustness studies 


8.5. Conclusion on robustness studies 


A robustness evaluation strategy was developed, which uses the results of the 
optimization as input data. Therefore, a data selection procedure has been 
presented, which focuses on a selection of the training data considering that 
the meta-model has to be more precise around the Pareto front, since the 
robustness is evaluated for those points. 


The robustness was evaluated taking two different approaches into account. 
First, the degree of robustness calculates the maximum possible changes of 
the design parameters, which result in solutions that are still robust. The sec- 
ond approach is the robustness index, where the worst case deviation of the 
design parameters is evaluated to check whether it is still robust or not. 


The approach is demonstrated with a heat map created for this purpose. The 
heat map is a useful tool for process and quality engineers, since it shows the 
areas in which the patch has to be, to remain the optimal solution. This 
means, the heat map can be used as an additional decision making tool. 


133 


9 Conclusion and outlook 


9.1. Conclusion 


In this work an optimization strategy for the positioning of local fiber rein- 
forcement patches, under consideration of manufacturing constraints, has 
been developed. Therefore an evolutionary algorithm has been used as opti- 
mization algorithm. The principal concept, consisting of the steps fitness eval- 
uation, crossover, mutation, and selection has been discussed, and their in- 
fluence on the optimization performance has been studied. To consider 
manufacturing effects in conjunction with the structural performance, a kin- 
ematic draping, a warpage, and a structural simulation are included in the 
fitness evaluation of the workflow. 


The first step is to perform a kinematic draping simulation to predict the final, 
deformed geometry of the patch. A kinematic draping simulation has been 
chosen, since a large number of function evaluations is necessary during the 
optimization process, and therefore a very efficient method is required. 
Therefore, a comparison of kinematic and mechanical draping methods has 
been presented to show the advantages and disadvantages of a kinematic 
draping approach. Within this study, three patch configurations with increas- 
ing degree of geometrical complexity have been compared on the basis of a 
three-dimensionally curved reference structure. Taking the results of the 
comparison with the mechanical draping simulation into account, the kine- 
matic draping simulation can be seen as a good approximation of the patch’s 
forming behavior. Compared with the mechanical method, the kinematic 
draping simulation is much faster, and requires no material parameters. 
Therefore, the kinematic draping suits best for the proposed optimization 
workflow. Additionally, the optimization shall work as a suggestion for the 
designer and not as a final verification. Therefore, a mechanical draping sim- 
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ulation must be performed afterwards to ensure manufacturability of the op- 
timized patch setup without forming defects like wrinkling, fiber fracture or 
gapping. 


For the warpage simulation, the curing behavior of the resin has to be pre- 
dicted. Therefore, an approach presented by Kamal [73] was used. A compar- 
ison of the simulation results with experimental results, for the cure rate and 
degree of cure, showed a good accordance by different heat rates. Further- 
more, a process-induced strain calculation was performed based on the 
method proposed by Svanberg and Holmberg [82]. With their approach, the 
influence of the incorporation of chemical strain during the curing process 
could be demonstrated. Furthermore, the effective coefficients of thermal ex- 
pansion for the UD reinforcement patch were calculated in dependence on 
fiber orientation and fiber volume content with the method of Schapery 
[122]. The anisotropic thermal coefficients for the discontinuous fiber rein- 
forced component were calculated with the approach proposed by Rosen and 
Hashin [123]. Furthermore, the assumption was made, that the effective co- 
efficients of chemical shrinkage can be calculated with the same approach as 
the coefficients of thermal expansion. 


Beside the effects resulting from warpage and draping, further manufacturing 
inaccuracies influence the performance of the final part. Therefore, two main 
sources of deviations, orientation and size changes of the patch, could be 
identified and their consideration during the optimization process has been 
discussed. Thereby, orientation changes refer to misalignments of the whole 
patch or the fibers, whereas size changes occur in terms of length and width 
changes. Furthermore, the most important question to be answered is, how 
accurate the process has to be, to ensure that a predicted optimum remains 
an optimum for the real part. Therefore, a robustness evaluation is performed 
after the optimization. 


To rate the structural performance of a configuration a structural simulation 
is performed, in which the results from the draping simulation are included. 
Therefore, a linear elastic simulation of the component is performed, and the 
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strain energy is used to rate the fitness of the current configuration. For the 
evaluation of the optimizations performance and results, a number of criteria 
have been developed. The presented visualization approach by means of a 
heat map, is a helpful tool for the designer to decide, in which part of the 
component the reinforcement has the highest impact. 


The presented workflow and the evaluation criteria have been applied to a 
number of demonstration applications. Thereby, the degree of complexity 
has been increased from a simple plate structure to a complex 3D component. 
First, the plate structure was used to perform basic feasibility studies and to 
compare the performance of the different crossover operators. Thereby, the 
multi-point binary crossover offers a better performance in terms of conver- 
gence rate and number of found solutions, while the found Pareto fronts are 
similar with both crossover types. Utilizing an optimization run with a large 
number of generations, the ability of the dominated area to work as a con- 
vergence criteria has been demonstrated. While increasing the complexity of 
the application example, the influence of different settings of the evolution- 
ary algorithm on the optimization’s performance has been studied. Thereby, 
all different settings converge to a similar Pareto front, which shows that the 
approach is capable to trustworthy solve the patch optimization problem. 
Thereby, the use of a niching pressure ensures a more balanced result, while 
a higher mutation rate results in a better coverage of the extreme areas of 
the search space. The last application example studied here, was the refer- 
ence structure developed within the DFG-GRK2078. This example showed the 
ability of the proposed approach to also solve complex optimization prob- 
lems. Furthermore, the reference structure has been used to perform an op- 
timization with the three objectives, strain energy from a structural load case, 
and distortion and strain energy from curing. Thereby, could be shown that it 
is possible to find a three-dimensional Pareto front with the proposed ap- 
proach. Furthermore, it was shown that the combination of distortion and 
strain energy from curing, is a suitable approach to consider manufacturing 
influences in the optimization. 
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Based on the final Pareto front a decision for one patch configuration has to 
be made by the designer. Since the presented optimization workflow results 
in a number of Pareto optimal solutions, usually a trade-off between the ob- 
jective functions has to be made when choosing one solution. Therefore, the 
robustness of a solution is introduced as an additional decision criterion. The 
developed robustness study uses the results of the optimization as input data, 
to create a meta model, which is then used during the robustness evaluation. 
Therefore, a data selection procedure has been presented, which focuses on 
a selection of the training data considering that the meta-model has to be 
more precise around the Pareto front. The robustness was evaluated taking 
two different approaches into account. First, the degree of robustness calcu- 
lates the maximum possible changes of the design parameters, which result 
in solutions that are still robust. The second approach is the robustness index, 
where the worst case deviation of the design parameters is evaluated to 
check whether it is still robust or not. Finally, the proposed robustness work- 
flow was demonstrated, using the slightly curved structure. Thereby could be 
shown that the robustness is a helpful further criteria for the decision maker 
for the selection of a solution on the final Pareto front. Additionally, the ro- 
bustness measurements are useful tool for the consideration of manufactur- 
ing inaccuracies during the part design process. 


9.2. Outlook 


The proposed optimization workflow showed its ability to optimize local rein- 
forcement patches. Furthermore, the integration of manufacturing con- 
straints as well as draping and warpage simulation within the optimization 
workflow improved the quality of the results. Here, an integration of the mold 
filling process in the optimization loop could be useful, to further improve the 
quality of the results. Thereby, the influence of the mold filling process on the 
final patch position and geometry could be considered during the optimiza- 
tion. 
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The robustness study showed that the use of meta models, instead of com- 
putationally expensive simulations, is a promising approach. Hence, the opti- 
mization performance could be enhanced by integrating meta models in the 
optimization workflow. Therefore, two different approaches seem to be 
promising, the setup of meta models a-priori or an alternating use of meta 
models and simulations. By using a-priori created meta models, a number of 
simulations is performed, according to a design of experiments (DoE), to cre- 
ate the required training data. For the second approach, meta models and 
simulation are used alternating, where the data created by the simulations 
can be used to train the meta models. Thereby, the meta models could espe- 
cially be used to replace the warpage simulation, since it is computationally 
more expensive than the structural simulation. 


In addition, the robustness evaluation could be directly integrated into the 
optimization loop. Thereby, the robustness of a solution could serve either as 
constraint or as and additional objective of the optimization. 
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Diese Arbeit leistet einen Beitrag zur Optimierung von lokalen endlosfaserver- 
starkten Patches, unter Berücksichtigung von Fertigungsbedingungen. 

Der in dieser Arbeit entwickelte Mehrziel-Optimierungsansatz verwendet dabei 
einen evolutionaren Algorithmus als Optimierungsalgorithmus. Da diese Klasse 
von Algorithmen eine Vielzahl von Funktionsauswertungen erfordert, sind effi- 
ziente Methoden zur Bewertung der Fitnesswerte notwendig. Aus diesem 
Grund wird eine kinematische Drapiersimulation zur Vorhersage der Patch- 
Geometrie verwendet. Als Zielfunktionen für die Mehrzieloptimierung werden 
sowohl strukturelle als auch prozessbezogene Ziele verwendet. Als strukturelles 
Optimierungsziel wird die Bauteilsteifigkeit verwendet, welche mittels einer li- 
near-elastischen Struktursimulation ermittelt wird, während für das prozessbe- 
zogene Ziel eine Verzugssimulation durchgeführt wird. Neben dem Verzugsziel 
werden weitere relevante Fertigungseinflüsse und -randbedingungen aus dem 
Halbzeug-, Drapier- und Co-Molding-Prozess diskutiert und finden während 
der Optimierung Berücksichtigung. 
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